Mercurial > repos > public > sbplib
annotate +time/Rungekutta4SecondOrder.m @ 774:66eb4a2bbb72 feature/grids
Remove default scaling of the system.
The scaling doens't seem to help actual solutions. One example that fails in the flexural code.
With large timesteps the solutions seems to blow up. One particular example is profilePresentation
on the tdb_presentation_figures branch with k = 0.0005
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Wed, 18 Jul 2018 15:42:52 -0700 |
parents | ae905a11e32c |
children | c6fcee3fcf1b 8894e9c49e40 74eec7e69b63 |
rev | line source |
---|---|
0 | 1 classdef Rungekutta4SecondOrder < 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 | |
549
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
18 % Solves u_tt = Du + Eu_t + S by |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
19 % Rewriting on first order form: |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
20 % w_t = M*w + C(t) |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
21 % where |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
22 % M = [ |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
23 % 0, I; |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
24 % D, E; |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
25 % ] |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
26 % and |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
27 % C(t) = [ |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
28 % 0; |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
29 % S(t) |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
30 % ] |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
31 % D, E, S can either all be constants or all be function handles, |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
32 % They can also be omitted by setting them equal to the empty matrix. |
0 | 33 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) |
34 obj.D = D; | |
35 obj.E = E; | |
36 obj.S = S; | |
37 obj.m = length(v0); | |
13
b18d3d201a71
Fixed initialization of step counter in timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
38 obj.n = 0; |
0 | 39 |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
40 |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
41 if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle') |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
42 default_arg('D', @(t)sparse(obj.m, obj.m)); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
43 default_arg('E', @(t)sparse(obj.m, obj.m)); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
44 default_arg('S', @(t)sparse(obj.m, 1) ); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
45 |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
46 if ~isa(D, 'function_handle') |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
47 D = @(t)D; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
48 end |
345
7b5ef8b89268
Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents:
222
diff
changeset
|
49 if ~isa(E, 'function_handle') |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
50 E = @(t)E; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
51 end |
345
7b5ef8b89268
Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents:
222
diff
changeset
|
52 if ~isa(S, 'function_handle') |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
53 S = @(t)S; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
54 end |
0 | 55 |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
56 obj.k = k; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
57 obj.t = t0; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
58 obj.w = [v0; v0t]; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
59 |
549
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
60 % Avoid matrix formulation because it is VERY slow |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
61 obj.F = @(w,t)[ |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
62 w(obj.m+1:end); |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
63 D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t); |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
64 ]; |
0 | 65 else |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
66 |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
67 default_arg('D', sparse(obj.m, obj.m)); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
68 default_arg('E', sparse(obj.m, obj.m)); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
69 default_arg('S', sparse(obj.m, 1) ); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
70 |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
71 I = speye(obj.m); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
72 O = sparse(obj.m,obj.m); |
0 | 73 |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
74 obj.M = [ |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
75 O, I; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
76 D, E; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
77 ]; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
78 obj.C = [ |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
79 zeros(obj.m,1); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
80 S; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
81 ]; |
0 | 82 |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
83 obj.k = k; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
84 obj.t = t0; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
85 obj.w = [v0; v0t]; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
86 |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
87 obj.F = @(w,t)(obj.M*w + obj.C); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
88 end |
0 | 89 end |
90 | |
91 function [v,t] = getV(obj) | |
92 v = obj.w(1:end/2); | |
93 t = obj.t; | |
94 end | |
95 | |
96 function [vt,t] = getVt(obj) | |
97 vt = obj.w(end/2+1:end); | |
98 t = obj.t; | |
99 end | |
100 | |
101 function obj = step(obj) | |
102 obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); | |
103 obj.t = obj.t + obj.k; | |
104 obj.n = obj.n + 1; | |
105 end | |
106 end | |
107 | |
108 | |
109 methods (Static) | |
110 function k = getTimeStep(lambda) | |
111 k = rk4.get_rk4_time_step(lambda); | |
112 end | |
113 end | |
114 | |
115 end |