Mercurial > repos > public > sbplib
changeset 549:ae905a11e32c feature/grids
Change away from matrix formulation in Rk4 second order form
The new formulation can be up to 8 times faster
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Mon, 28 Aug 2017 11:06:00 +0200 |
| parents | 24b2487b01c2 |
| children | e860670e72f1 |
| files | +time/Rungekutta4SecondOrder.m |
| diffstat | 1 files changed, 18 insertions(+), 14 deletions(-) [+] |
line wrap: on
line diff
diff -r 24b2487b01c2 -r ae905a11e32c +time/Rungekutta4SecondOrder.m --- a/+time/Rungekutta4SecondOrder.m Thu Aug 24 11:05:31 2017 +0200 +++ b/+time/Rungekutta4SecondOrder.m Mon Aug 28 11:06:00 2017 +0200 @@ -15,7 +15,19 @@ methods - % Solves u_tt = Du + Eu_t + S + % Solves u_tt = Du + Eu_t + S by + % Rewriting on first order form: + % w_t = M*w + C(t) + % where + % M = [ + % 0, I; + % D, E; + % ] + % and + % C(t) = [ + % 0; + % S(t) + % ] % D, E, S can either all be constants or all be function handles, % They can also be omitted by setting them equal to the empty matrix. function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) @@ -41,23 +53,15 @@ S = @(t)S; end - I = speye(obj.m); - O = sparse(obj.m,obj.m); - - obj.M = @(t)[ - O, I; - D(t), E(t); - ]; - obj.C = @(t)[ - zeros(obj.m,1); - S(t); - ]; - obj.k = k; obj.t = t0; obj.w = [v0; v0t]; - obj.F = @(w,t)(obj.M(t)*w + obj.C(t)); + % Avoid matrix formulation because it is VERY slow + obj.F = @(w,t)[ + w(obj.m+1:end); + D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t); + ]; else default_arg('D', sparse(obj.m, obj.m));