|
0
|
1 classdef Ode45 < time.Timestepper
|
|
|
2 properties
|
|
|
3 F
|
|
|
4 k
|
|
|
5 t
|
|
|
6 w
|
|
|
7 m
|
|
|
8 D
|
|
|
9 E
|
|
|
10 S
|
|
|
11 M
|
|
|
12 C
|
|
|
13 n
|
|
|
14 end
|
|
|
15
|
|
|
16
|
|
|
17 methods
|
|
|
18 function obj = Ode45(D, E, S, k, t0, v0, v0t)
|
|
|
19 obj.D = D;
|
|
|
20 obj.E = E;
|
|
|
21 obj.S = S;
|
|
|
22 obj.m = length(v0);
|
|
|
23
|
|
|
24 I = speye(obj.m);
|
|
|
25 O = sparse(obj.m,obj.m);
|
|
|
26 obj.M = [O, I; D, E*I]; % Multiply with I to allow 0 as input.
|
|
|
27
|
|
|
28 if S == 0
|
|
|
29 obj.C = zeros(2*obj.m,1);
|
|
|
30 else
|
|
|
31 obj.C = [zeros(obj.m,1), S];
|
|
|
32 end
|
|
|
33
|
|
|
34 obj.k = k;
|
|
|
35 obj.t = t0;
|
|
|
36 obj.w = [v0; v0t];
|
|
|
37
|
|
|
38 obj.F = @(w,t)(obj.M*w + obj.C);
|
|
|
39 end
|
|
|
40
|
|
|
41 function [v,t] = getV(obj)
|
|
|
42 v = obj.w(1:end/2);
|
|
|
43 t = obj.t;
|
|
|
44 end
|
|
|
45
|
|
|
46 function [vt,t] = getVt(obj)
|
|
|
47 vt = obj.w(end/2+1:end);
|
|
|
48 t = obj.t;
|
|
|
49 end
|
|
|
50
|
|
|
51 function obj = step(obj)
|
|
|
52 [t,w] = ode45(@(t,w)(obj.F(w,t)),[obj.t obj.t+obj.k],obj.w);
|
|
|
53
|
|
|
54 obj.t = t(end);
|
|
|
55 obj.w = w(end,:)';
|
|
|
56 obj.n = obj.n + 1;
|
|
|
57 end
|
|
|
58 end
|
|
|
59
|
|
|
60
|
|
|
61 methods (Static)
|
|
|
62 function k = getTimeStep(lambda)
|
|
|
63 k = rk4.get_rk4_time_step(lambda);
|
|
|
64 end
|
|
|
65 end
|
|
|
66
|
|
|
67 end |
