Mercurial > repos > public > sbplib

comparison +time/CdiffNonlin.m @ 1:5ae4f23d9130

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 17 Sep 2015 19:19:58 +0200
parents 48b6fb693025
children 8e14b5a577a6
comparison
equal deleted inserted replaced
0:48b6fb693025 1:5ae4f23d9130
1 classdef Cdiff < time.Timestepper 1 classdef CdiffNonlin < time.Timestepper
2 properties 2 properties
3 D 3 D
4 E 4 E
5 S 5 S
6 k 6 k
10 n 10 n
11 end 11 end
12 12
13 13
14 methods 14 methods
15 function obj = Cdiff(D, E, S, k, t0, v, v_prev) 15 function obj = CdiffNonlin(D, E, S, k, t0, v, v_prev)
16 m = size(D,1); 16 default_arg('S',0);
17 default_arg('E',sparse(m,m)); 17 default_arg('E',0);
18 default_arg('S',sparse(m,1));
19 18
20 if ~(issparse(D) && issparse(E) && issparse(S)) 19 if isnumeric(S) && S == 0
21 warning('One of the matrices D, E, S is not sparse!') 20 S = @(v)0;
22 print_issparse(D)
23 print_issparse(E)
24 print_issparse(S)
25 end 21 end
22
23 if isnumeric(E) && E == 0
24 E = @(v)0;
25 end
26
27
28 % m = size(D,1);
29 % default_arg('E',sparse(m,m));
30 % default_arg('S',sparse(m,1));
26 31
27 obj.D = D; 32 obj.D = D;
28 obj.E = E; 33 obj.E = E;
29 obj.S = S; 34 obj.S = S;
30 obj.k = k; 35 obj.k = k;
31 obj.t = t0+k; 36 obj.t = t0;
32 obj.v = v; 37 obj.v = v;
33 obj.v_prev = v_prev; 38 obj.v_prev = v_prev;
34 end 39 end
35 40
36 function [v,t] = getV(obj) 41 function [v,t] = getV(obj)
42 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) 47 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
43 t = obj.t; 48 t = obj.t;
44 end 49 end
45 50
46 function obj = step(obj) 51 function obj = step(obj)
47 [obj.v, obj.v_prev] = time.cdiff.cdiff(obj.v, obj.v_prev, obj.k, obj.D, obj.E, obj.S); 52 [obj.v, obj.v_prev] = time.cdiff.cdiff(obj.v, obj.v_prev, obj.k, obj.D(obj.v), obj.E(obj.v), obj.S(obj.v));
48 obj.t = obj.t + obj.k; 53 obj.t = obj.t + obj.k;
49 obj.n = obj.n + 1; 54 obj.n = obj.n + 1;
50 end 55 end
51 end 56 end
52 end 57 end