Mercurial > repos > public > sbplib

comparison +time/CdiffNonlin.m @ 31:d1f9dd55a2b0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed bug in step() for CdiffNonlin.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 28 Sep 2015 21:56:02 +0200
parents cd2e28c5ecd2
children b5e5b195da1e
comparison
equal deleted inserted replaced
30:cd2e28c5ecd2 31:d1f9dd55a2b0
65 65
66 66
67 %% Calculate matrices need for the timestep 67 %% Calculate matrices need for the timestep
68 % Before optimization: A = 1/k^2 * I - 1/(2*k)*E; 68 % Before optimization: A = 1/k^2 * I - 1/(2*k)*E;
69 k = obj.k; 69 k = obj.k;
70
70 Aj = 1/k^2 * I(j,j) - 1/(2*k)*E(j,j); 71 Aj = 1/k^2 * I(j,j) - 1/(2*k)*E(j,j);
71 B = 2/k^2 * I + D; 72 B = 2/k^2 * I + D;
72 C = -1/k^2 * I - 1/(2*k)*E; 73 C = -1/k^2 * I - 1/(2*k)*E;
73 74
74 %% Take the timestep 75 %% Take the timestep
75 v = obj.v; 76 v = obj.v;
77 v_prev = obj.v_prev;
76 78
77 % Before optimization: obj.v = A\(B*v + C*v_prev + S); 79 % Want to solve the seq A*v_next = b where
78 obj.v(i) = k^2*(B(i,i)*v(i) -1/k^2*obj.v_prev(i) + S(i)); 80 b = (B*v + C*v_prev + S);
79 obj.v(j) = Aj\(B(j,j)*v(j) + C(j,j)*obj.v_prev(j) + S(j)); 81
82 % Before optimization: obj.v = A\b;
83
84 obj.v(i) = k^2*b(i);
85 obj.v(j) = Aj\b(j);
86
80 obj.v_prev = v; 87 obj.v_prev = v;
81 88
82 %% Update state of the timestepper 89 %% Update state of the timestepper
83 obj.t = obj.t + obj.k; 90 obj.t = obj.t + obj.k;
84 obj.n = obj.n + 1; 91 obj.n = obj.n + 1;