Mercurial > repos > public > sbplib

comparison +time/Ode45.m @ 0:48b6fb693025

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial commit.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 17 Sep 2015 10:12:50 +0200
parents
children b5e5b195da1e
comparison
equal deleted inserted replaced
-1:000000000000 0:48b6fb693025
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