Problem with finite difference
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Paul Eluard
am 12 Nov. 2017
Beantwortet: David Goodmanson
am 13 Nov. 2017
Hello !
In order to solve the following ODE : $U_{xx} + U_x - 2U + 2 = 0$ with boundary condition $U(-1) = 0 = U(1)$, I use centered differences of order 2 :
$U_{xx} = \frac{U_{i+1} - 2 U_i + U_{i-1}}{h^2} + O(h^2)$
$U_x = \frac{U_{i+1} - U_{i-1}}{2h} + O(h^2)$
Putting this altogether gives me the following code :
Mxx = gallery('tridiag',n+1,1,-2,1)*1/h^2; Mx = gallery('tridiag',n+1,1,0,-1)*1/(2*h); M0 = speye(n+1); M = Mxx + Mx - 2*M0; M(1,1) = 1; M(1,2) = 0; M(end,end) = 1; M(end,end-1) = 0; M = sparse(M); B = [0;-2*ones(n-1,1);0]; U = M\B;
But the solution doesn't fit the analytical one at all... $u(x) = 1 - \frac{\sinh(2)e^x + \sinh(1)e^{-2x}}{\sinh(3)}$
![](/matlabcentral/answers/uploaded_files/94340/untitled.jpg)
Where is the mistake ?
0 Kommentare
Akzeptierte Antwort
David Goodmanson
am 13 Nov. 2017
Hi Paul,
You did not say what h is, but I took it to be 2/n. After replacing Mx with its transpose, which changes the sign of the derivative (the results say it's now the correct sign but I didn't sit down and prove it) you get really good agreement. The plot is basically an overlay.
n = 1000;
h=2/n;
Mxx = gallery('tridiag',n+1,1,-2,1)*1/h^2;
Mx = gallery('tridiag',n+1,1,0,-1)'*1/(2*h); % transposed
M0 = speye(n+1);
M = Mxx + Mx - 2*M0;
M(1,1) = 1; M(1,2) = 0;
M(end,end) = 1; M(end,end-1) = 0;
M = sparse(M);
B = [0;-2*ones(n-1,1);0];
U = M\B;
x = linspace(-1,1,1001);
u = 1 - (sinh(2)*exp(x) + sinh(1)*exp(-2*x))/sinh(3);
figure(1)
plot(1:n+1,U',1:n+1,u)
max(abs(U'-u))
0 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!