-
Notifications
You must be signed in to change notification settings - Fork 0
/
RK.m
99 lines (95 loc) · 2.93 KB
/
RK.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
function [t,x,y] = RK(method,dx,dy,xi,yi,ti,tf,h)
% Runge-Kutta method for solving the 2nd order ODE y'' = f(t,x,y)
% This is achieved by writing the 2nd order ODE as two coupled 1st order ODEs:
% dx = y = x'
% dy = y' = x''
% then solving both together using a specified RK method.
%
% Parameters:
% method: name of the RK method to be used (case-insensitive)
% dx, dy: coupled 1st order ODEs equivalent of the 2nd order ODE y'' = f(t,x,y)
% dx = y = x'
% dy = y' = x''
% xi: initial value of x
% yi: initial value of y
% ti: initial value of t
% tf: final value of t
% h: time interval size
% Generate the necessary matrices:
% s: number of stages = order of the method
% a: Runge-Kutta matrix
% b: weights
% c: nodes
if(strcmpi(method,'euler')) % Euler's method, RK1
s = 1;
a = 0;
b = 1;
c = 0;
elseif(strcmpi(method,'midpoint')) % Explicit midpoint method, RK2
s = 2;
a = [0 0;
1/2 0];
b = [0 1];
c = [0 1/2];
elseif(strcmpi(method,'heun')) % Heun's method, RK2
s = 2;
a = [0 0;
1 0];
b = [1/2 1/2];
c = [0 1];
elseif(strcmpi(method,'ralston')) % Ralston's method, RK2
s = 2;
a = [0 0;
2/3 0];
b = [1/4 3/4];
c = [0 2/3];
elseif(strcmpi(method,'rk3')) % Kutta's 3rd order method, RK3
s = 3;
a = [0 0 0;
1/2 0 0;
-1 2 0];
b = [1/6 2/3 1/6];
c = [0 1/2 1];
elseif(strcmpi(method,'rk4')) % Classic 4th order method, RK4
s = 4;
a = [0 0 0 0;
1/2 0 0 0;
0 1/2 0 0;
0 0 1 0];
b = [1/6 1/3 1/3 1/6];
c = [0 1/2 1/2 1];
elseif(strcmpi(method,'3/8')) % 3/8 4th order method, RK4
s = 4;
a = [0 0 0 0;
1/3 0 0 0;
-1/3 1 0 0;
1 -1 1 0];
b = [1/8 3/8 3/8 3/8];
c = [0 1/3 2/3 1];
else
error('Method cannot be found.');
end
% Applying numerical method
N = round((tf-ti)/h); % total number of steps
t = zeros(1,N); % [1 x N] matrix of 0's
x = zeros(1,N);
y = zeros(1,N);
% Set initial conditions
t(1) = ti;
x(1) = xi;
y(1) = yi;
% Iteration
for n = 1:N-1
k = zeros(1,s);
l = zeros(1,s);
for j = 1:s
tt = t(n) + h * c(j);
xx = x(n) + dot(a(j,:),k);
yy = y(n) + dot(a(j,:),l);
k(j) = h * dx(tt, xx, yy);
l(j) = h * dy(tt, xx, yy);
end
t(n+1) = t(n) + h;
x(n+1) = x(n) + dot(b,k);
y(n+1) = y(n) + dot(b,l);
end