-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathrk4.h
117 lines (105 loc) · 2.65 KB
/
rk4.h
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
/*
* rk4.h
* Copyright (C) Farooq Mela 2003-2013.
*/
#ifndef _RK4_H_
#define _RK4_H_
/* Fourth-order Runge-Kutta step, using array of functions.
*
* Input parameters:
* n: number of equations
* f: array of pointers to differential equations
* t: current time parameter
* h: time step or delta
* y: value of y at time t
* Output parameter:
* y_next: value of y at time t + h
*/
void
rk4(const int n, double (*f[])(double t, const double y[]),
double t, double h, const double y[], double y_next[])
{
/* k1 = h * f(t, y) */
double k[n];
for (int i = 0; i < n; ++i) {
k[i] = h * (*f[i])(t, y);
y_next[i] = k[i];
}
/* k2 = h * f(t + h/2, y + k1/2) */
double temp[n];
for (int i = 0; i < n; ++i)
temp[i] = y[i] + k[i]/2;
for (int i = 0; i < n; ++i) {
k[i] = h * (*f[i])(t + h/2, temp);
y_next[i] += k[i] * 2;
}
/* k3 = h * f(t + h/2, y + k2/2) */
for (int i = 0; i < n; ++i)
temp[i] = y[i] + k[i]/2;
for (int i = 0; i < n; ++i) {
k[i] = h * (*f[i])(t + h/2, temp);
y_next[i] += k[i] * 2;
}
/* k4 = h * f(t + h, y + k3) */
for (int i = 0; i < n; ++i)
temp[i] = y[i] + k[i];
for (int i = 0; i < n; ++i) {
k[i] = h * (*f[i])(t + h, temp);
y_next[i] += k[i];
}
/* y_next = y + (k1 + 2(k2 + k3) + k4) / 6 */
for (int i = 0; i < n; ++i)
y_next[i] = y[i] + y_next[i] / 6.;
}
/* Fourth-order Runge-Kutta step, using vector-valued function.
*
* Input parameters:
* n: number of equations
* f: pointer to vector-valued differential equation
* t: current time parameter
* h: time step or delta
* y: value of y at time t
* Output parameter:
* y_next: value of y at time t + h
*/
void
rk4v(const int n, void (*f)(double t, const double y[], double dy[]),
double t, double h, const double y[], double y_next[])
{
/* k1 = h * f(t, y) */
double k[n];
f(t, y, k);
for (int i = 0; i < n; ++i) {
k[i] *= h;
y_next[i] = k[i];
}
/* k2 = h * (t + h/2, y + k1/2) */
double temp[n];
for (int i = 0; i < n; ++i)
temp[i] = y[i] + k[i]/2;
f(t + h/2, temp, k);
for (int i = 0; i < n; ++i) {
k[i] *= h;
y_next[i] += k[i] * 2;
}
/* k3 = h * (t + h/2, y + k2/2) */
for (int i = 0; i < n; ++i)
temp[i] = y[i] + k[i]/2;
f(t + h/2, temp, k);
for (int i = 0; i < n; ++i) {
k[i] *= h;
y_next[i] += k[i] * 2;
}
/* k4 = h * (t + h, y + k3) */
for (int i = 0; i < n; ++i)
temp[i] = y[i] + k[i];
f(t + h, temp, k);
for (int i = 0; i < n; ++i) {
k[i] *= h;
y_next[i] += k[i];
}
/* y_next = y + (k1 + 2(k2 + k3) + k4) / 6 */
for (int i = 0; i < n; ++i)
y_next[i] = y[i] + y_next[i] / 6.;
}
#endif