1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
#include "diffeq.h"
#include <string.h>
#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
/* Calling:
* if data is non-null, t is used. if data is null, t is not used.
*
* 1. Initializing files: files = [array of files], data = NULL, n = sizeof [array of files]
* 2. Updating files: files = anything non-NULL [unused], data = [array of data], n = sizeof [array of data].
* n may be 0, then saved n from init is used
* 3. Closing files: files = NULL, data = NULL, n = anything
*
* returns: whether the operation was successful
*/
char flush_approx(const char ** files, f64 t, f64 * data, usz n) {
static FILE ** handlers;
static usz size;
usz i;
if (data == NULL && files == NULL) {
for (i = 0; i < n; ++i) {
if (fclose(handlers[i]) != 0)
perror("closing files");
}
free(handlers);
return 1;
}
if (data == NULL) {
size = n;
handlers = calloc(n, sizeof (FILE*));
for (i = 0; i < n; ++i) {
handlers[i] = fopen(files[i], "w+");
if (handlers[i] == NULL) {
perror(files[i]);
return 0;
}
}
return 1;
}
if (files == NULL) {
for (i = 0; i < (n == 0 ? size : n); ++i)
fprintf(handlers[i], "%f %f\n", t, data[i]);
return 1;
}
return 0;
}
void difeq_solve_euler(eqf_t *eqfs, f64 *init, const char **files, usz n, f64 a, f64 b, f64 dt) {
f64 *guesses1 = calloc(n, sizeof(f64));
f64 *guesses2 = calloc(n, sizeof(f64));
f64 cur = a;
size_t i;
flush_approx(files, 0, NULL, n);
memcpy(guesses1, init, n * sizeof(f64));
memcpy(guesses2, init, n * sizeof(f64));
for(cur = a; cur < b; cur += dt) {
for (i = 0; i < n; ++i)
guesses2[i] += dt * eqfs[i](guesses1, n, cur);
memcpy(guesses1, guesses2, n * sizeof(f64));
flush_approx(NULL, cur, guesses1, n);
}
flush_approx(NULL, 0, NULL, 0);
free(guesses1);
free(guesses2);
}
void difeq_solve_RK(eqf_t *eqfs, f64 *init, const char **files, usz n, f64 a, f64 b, f64 dt) {
f64 *guesses1 = calloc(n, sizeof(f64));
f64 *guesses2 = calloc(n, sizeof(f64));
f64 cur = a;
f64 *k1s = calloc(n, sizeof(f64)),
*k2s = calloc(n, sizeof(f64)),
*k3s = calloc(n, sizeof(f64)),
*k4s = calloc(n, sizeof(f64));
size_t step, i;
flush_approx(files, 0, NULL, n);
memcpy(guesses2, init, n * sizeof(f64));
memcpy(guesses1, init, n * sizeof(f64));
for(step = 0; step < (b-a)/dt; cur = a + ++step*dt) {
for (i = 0; i < n; ++i)
k1s[i] = eqfs[i](guesses1, n, cur);
for (i = 0; i < n; ++i)
guesses2[i] = guesses1[i] + dt/2 * k1s[i];
for (i = 0; i < n; ++i)
k2s[i] = eqfs[i](guesses2, n, cur + dt/2);
for (i = 0; i < n; ++i)
guesses2[i] = guesses1[i] + dt/2 * k2s[i];
for (i = 0; i < n; ++i)
k3s[i] = eqfs[i](guesses2, n, cur + dt/2);
for (i = 0; i < n; ++i)
guesses2[i] = guesses1[i] + dt * k3s[i];
for (i = 0; i < n; ++i)
k4s[i] = eqfs[i](guesses2, n, cur + dt);
for (i = 0; i < n; ++i)
guesses2[i] += dt/6 * (k1s[i] + 2*k2s[i] +2*k3s[i] + k4s[i]);
memcpy(guesses1, guesses2, n * sizeof(f64));
flush_approx(NULL, cur, guesses1, n);
}
flush_approx(NULL, 0, NULL, 0);
free(guesses1);
free(guesses2);
free(k1s);
free(k2s);
free(k3s);
free(k4s);
}
|
