Mercurial > repos > public > sbplib

annotate +time/+rk/ExplicitSecondOrder.m @ 1102:d4c895d4b524 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add skeleton for time.rk.ExplicitSecondOrder
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 09 Apr 2019 22:17:07 +0200
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1102
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 classdef ExplicitSecondOrder < time.Timestepper
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 F % RHS of the ODE
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 dt % Time step
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 t % Time point
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 v, vt % Solution state
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 n % Time level
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 bt
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 methods
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 % Timesteps v_tt = F(t,v,vt), using the specified ButcherTableau
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 % from t = t0 with timestep dt and initial conditions v(0) = v0
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 function obj = ExplicitSecondOrder(F, dt, t0, v0, v0t, bt)
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 assertType(bt, 'time.rk.ButcherTableau')
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17 obj.F = F;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 obj.dt = dt;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 obj.t = t0;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 obj.v = v0;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21 obj.vt = v0t;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 obj.n = 0;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 assert(bt.isExplicit())
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25 obj.bt = bt;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 function [v,t] = getV(obj)
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29 v = obj.v;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30 t = obj.t;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 function [vt,t] = getVt(obj)
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 vt = obj.vt;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 t = obj.t;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
38 function obj = step(obj)
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39 s = obj.bt.nStages();
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40 a = obj.bt.a;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41 b = obj.bt.b;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42 c = obj.bt.c;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44 t = obj.t;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 v = obj.v;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46 vt = obj.vt;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47 dt = obj.dt;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 k1 = obj.F(t, v, v_t);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50 k2 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t, v_t + 1/2*dt*k1);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 k3 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t + 1/4*dt^2*k1, v_t + 1/2*dt*k2);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 k4 = obj.F(t + dt, v + dt*v_t + 1/2*dt^2*k2, v_t + dt*k3);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54 % Compute rates K
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55 K = zeros(length(v), s);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 for i = 1:s
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 U_i = obj.v;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 V_i = obj.vt;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 for j = 1:i-1
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60 U_i = U_i % + dt*a(i,j)*K(:,j);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
61 V_i = V_i % + dt*a(i,j)*K(:,j);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 K(:,i) = F(t+dt*c(i), U_i, V_i);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66 % Compute updated solution
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
67 v_next = v;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68 vt_next = vt;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69 for i = 1:s
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 v_next = v_next % + dt*b(i)*K(:,i);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 vt_next = vt_next % + dt*b(i)*K(:,i);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74 obj.v = v_next;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75 obj.vt = vt_next;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
76 obj.t = obj.t + obj.dt;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77 obj.n = obj.n + 1;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
81 % Returns a vector of time points, including substage points,
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82 % in the time interval [t0, tEnd].
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
84 function tvec = timePoints(obj, t0, tEnd)
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85 % TBD: Should this be implemented here or somewhere else?
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86 N = round( (tEnd-t0)/obj.dt );
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
87 tvec = zeros(N*obj.s, 1);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88 s = obj.coeffs.s;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89 c = obj.coeffs.c;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90 for i = 1:N
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 ind = (i-1)*s+1 : i*s;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92 tvec(ind) = ((i-1) + c')*obj.dt;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96 % Returns a vector of quadrature weights corresponding to grid points
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
97 % in time interval [t0, tEnd], substage points included.
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 function weights = quadWeights(obj, t0, tEnd)
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100 % TBD: Should this be implemented here or somewhere else?
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101 N = round( (tEnd-t0)/obj.dt );
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102 b = obj.coeffs.b;
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
103 weights = repmat(b', N, 1);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
104 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
105 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
106
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
107 methods(Static)
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
108 % TBD: Function name
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
109 function ts = methodFromStr(F, dt, t0, v0, methodStr)
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
110 try
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
111 bt = time.rk.ButcherTableau.(method);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
112 catch
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
113 error('Runge-Kutta method ''%s'' is not implemented', methodStr)
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
114 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
115
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
116 ts = time.rk.Explicit(F, dt, t0, v0, bt);
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
117 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
118 end
d4c895d4b524 Add skeleton for time.rk.ExplicitSecondOrder
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
119 end