-
Notifications
You must be signed in to change notification settings - Fork 0
/
rkf45.m
119 lines (106 loc) · 4.24 KB
/
rkf45.m
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
118
119
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function [tout, yout] = rkf45(ode_function, tspan, y0, tolerance)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%{
This function uses the Runge-Kutta-Fehlberg 4(5) algorithm to
integrate a system of first-order differential equations
dy/dt = f(t,y).
y - column vector of solutions
f - column vector of the derivatives dy/dt
t - time
a - Fehlberg coefficients for locating the six solution
points (nodes) within each time interval.
b - Fehlberg coupling coefficients for computing the
derivatives at each interior point
c4 - Fehlberg coefficients for the fourth-order solution
c5 - Fehlberg coefficients for the fifth-order solution
tol - allowable truncation error
ode_function - handle for user M-function in which the derivatives f
are computed
tspan - the vector [t0 tf] giving the time interval for the
solution
t0 - initial time
tf - final time
y0 - column vector of initial values of the vector y
tout - column vector of times at which y was evaluated
yout - a matrix, each row of which contains the components of y
evaluated at the correponding time in tout
h - time step
hmin - minimum allowable time step
ti - time at the beginning of a time step
yi - values of y at the beginning of a time step
t_inner - time within a given time step
y_inner - values of y witin a given time step
te - trucation error for each y at a given time step
te_allowed - allowable truncation error
te_max - maximum absolute value of the components of te
ymax - maximum absolute value of the components of y
tol - relative tolerance
delta - fractional change in step size
eps - unit roundoff error (the smallest number for which
1 + eps > 1)
eps(x) - the smallest number such that x + eps(x) = x
User M-function required: ode_function
%}
% ---------------------------------------------------------------
a = [0 1/4 3/8 12/13 1 1/2];
b = [ 0 0 0 0 0
1/4 0 0 0 0
3/32 9/32 0 0 0
1932/2197 -7200/2197 7296/2197 0 0
439/216 -8 3680/513 -845/4104 0
-8/27 2 -3544/2565 1859/4104 -11/40];
c4 = [25/216 0 1408/2565 2197/4104 -1/5 0 ];
c5 = [16/135 0 6656/12825 28561/56430 -9/50 2/55];
if nargin < 4
tol = 1.e-8;
else
tol = tolerance;
end
t0 = tspan(1);
tf = tspan(2);
t = t0;
y = y0;
tout = t;
yout = y';
h = (tf - t0)/100; % Assumed initial time step.
while t < tf
hmin = 16*eps(t);
ti = t;
yi = y;
%...Evaluate the time derivative(s) at six points within the current
% interval:
for i = 1:6
t_inner = ti + a(i)*h;
y_inner = yi;
for j = 1:i-1
y_inner = y_inner + h*b(i,j)*f(:,j);
end
f(:,i) = feval(ode_function, t_inner, y_inner);
end
%...Compute the maximum truncation error:
te = h*f*(c4' - c5'); % Difference between 4th and
% 5th order solutions
te_max = max(abs(te));
%...Compute the allowable truncation error:
ymax = max(abs(y));
te_allowed = tol*max(ymax,1.0);
%...Compute the fractional change in step size:
delta = (te_allowed/(te_max + eps))^(1/5);
%...If the truncation error is in bounds, then update the solution:
if te_max <= te_allowed
h = min(h, tf-t);
t = t + h;
y = yi + h*f*c5';
tout = [tout;t];
yout = [yout;y'];
end
%...Update the time step:
h = min(delta*h, 4*h);
if h < hmin
fprintf(['\n\n Warning: Step size fell below its minimum\n'...
' allowable value (%g) at time %g.\n\n'], hmin, t)
return
end
end
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~