Mercurial > repos > public > sbplib

annotate +time/CdiffImplicit.m @ 1223:9fddc8749445 rv_diffOp_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Closing branch
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 05 Aug 2019 10:48:37 +0200
parents d5bce13ece23
children d6ede7f5cbf9
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 classdef CdiffImplicit < time.Timestepper
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 A, B, C, G
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 AA, BB, CC
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 k
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 t
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 v, v_prev
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 n
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 % LU factorization
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 L,U,p,q
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 methods
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 % Solves
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 % A*u_tt + B*u + C*v_t = G(t)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17 % u(t0) = f1
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 % u_t(t0) = f2
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 % starting at time t0 with timestep k
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 function obj = CdiffImplicit(A, B, C, G, f1, f2, k, t0)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21 default_arg('A', []);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 default_arg('C', []);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 default_arg('G', []);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 default_arg('f1', 0);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25 default_arg('f2', 0);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 default_arg('t0', 0);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 m = size(B,1);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30 if isempty(A)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 A = speye(m);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 if isempty(C)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 C = sparse(m,m);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
38 if isempty(G)
380
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
39 G = @(t) sparse(m,1);
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42 if isempty(f1)
393
67ca8964f03c Fix index error in CdiffImplicit.
Jonatan Werpers <jonatan@werpers.com>
parents: 387
diff changeset
43 f1 = sparse(m,1);
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46 if isempty(f2)
393
67ca8964f03c Fix index error in CdiffImplicit.
Jonatan Werpers <jonatan@werpers.com>
parents: 387
diff changeset
47 f2 = sparse(m,1);
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50 obj.A = A;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 obj.B = B;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 obj.C = C;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53 obj.G = G;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55 AA = 1/k^2*A + 1/2*B + 1/(2*k)*C;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 BB = -2/k^2*A;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 CC = 1/k^2*A + 1/2*B - 1/(2*k)*C;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 % AA*v_next + BB*v + CC*v_prev == G(t_n)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60 obj.AA = AA;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
61 obj.BB = BB;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 obj.CC = CC;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 v_prev = f1;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65 I = speye(m);
454
d5bce13ece23 Switch to first order start in CdiffImplicit to make the scheme more robust.
Jonatan Werpers <jonatan@werpers.com>
parents: 393
diff changeset
66 % v = (1/k^2*A)\((1/k^2*A - 1/2*B)*f1 + (1/k*I - 1/2*C)*f2 + 1/2*G(0));
d5bce13ece23 Switch to first order start in CdiffImplicit to make the scheme more robust.
Jonatan Werpers <jonatan@werpers.com>
parents: 393
diff changeset
67 v = f1 + k*f2;
d5bce13ece23 Switch to first order start in CdiffImplicit to make the scheme more robust.
Jonatan Werpers <jonatan@werpers.com>
parents: 393
diff changeset
68
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 if ~issparse(A) || ~issparse(B) || ~issparse(C)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 error('LU factorization with full pivoting only works for sparse matrices.')
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74 [L,U,p,q] = lu(AA,'vector');
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
76 obj.L = L;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77 obj.U = U;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 obj.p = p;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79 obj.q = q;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
81
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82 obj.k = k;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83 obj.t = t0+k;
385
b361a04894e3 Fix bug in implicit cdiff. Add stepTo method for timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents: 383
diff changeset
84 obj.n = 1;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85 obj.v = v;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86 obj.v_prev = v_prev;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
87 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89 function [v,t] = getV(obj)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90 v = obj.v;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 t = obj.t;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94 function [vt,t] = getVt(obj)
387
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
95 % Calculate next time step to be able to do centered diff.
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
96 v_next = zeros(size(obj.v));
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
97 b = obj.G(obj.t) - obj.BB*obj.v - obj.CC*obj.v_prev;
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
98
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
99 y = obj.L\b(obj.p);
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
100 z = obj.U\y;
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
101 v_next(obj.q) = z;
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
102
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
103
45c69aff2f41 Swtich to using centerd difference for calculating the first derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 385
diff changeset
104 vt = (v_next-obj.v_prev)/(2*obj.k);
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
105 t = obj.t;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
106 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
107
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
108 function obj = step(obj)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
109 b = obj.G(obj.t) - obj.BB*obj.v - obj.CC*obj.v_prev;
380
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
110 obj.v_prev = obj.v;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
111
380
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
112 % % Backslash
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
113 % obj.v = obj.AA\b;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
114
380
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
115 % LU with column pivot
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
116 y = obj.L\b(obj.p);
383
151b08366d63 Fix bug in CdiffImplicit.
Jonatan Werpers <jonatan@werpers.com>
parents: 380
diff changeset
117 z = obj.U\y;
151b08366d63 Fix bug in CdiffImplicit.
Jonatan Werpers <jonatan@werpers.com>
parents: 380
diff changeset
118 obj.v(obj.q) = z;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
119
380
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
120 % Update time
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
121 obj.t = obj.t + obj.k;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
122 obj.n = obj.n + 1;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
123 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
124 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
125 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
126
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
127
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
128
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
129
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
130
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
131 %%% Derivation
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
132 % syms A B C G
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
133 % syms n k
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
134 % syms f1 f2
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
135
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
136 % v = symfun(sym('v(n)'),n);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
137
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
138
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
139 % d = A/k^2 * (v(n+1) - 2*v(n) +v(n-1)) + B/2*(v(n+1)+v(n-1)) + C/(2*k)*(v(n+1) - v(n-1)) == G
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
140 % ic1 = v(0) == f1
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
141 % ic2 = A/k*(v(1)-f1) + k/2*(B*f1 + C*f2 - G) - f2 == 0
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
142
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
143 % c = collect(d, [v(n) v(n-1) v(n+1)]) % (-(2*A)/k^2)*v(n) + (B/2 + A/k^2 - C/(2*k))*v(n - 1) + (B/2 + A/k^2 + C/(2*k))*v(n + 1) == G
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
144 % syms AA BB CC
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
145 % % AA = B/2 + A/k^2 + C/(2*k)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
146 % % BB = -(2*A)/k^2
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
147 % % CC = B/2 + A/k^2 - C/(2*k)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
148 % s = subs(c, [B/2 + A/k^2 + C/(2*k), -(2*A)/k^2, B/2 + A/k^2 - C/(2*k)], [AA, BB, CC])
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
149
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
150
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
151 % ic2_a = collect(ic2, [v(1) f1 f2]) % (A/k)*v(1) + ((B*k)/2 - A/k)*f1 + ((C*k)/2 - 1)*f2 - (G*k)/2 == 0
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
152