Mercurial > repos > public > sbplib

annotate +time/CdiffImplicit.m @ 774:66eb4a2bbb72 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 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