Problem with finite difference

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)}$

Where is the mistake ?

Answer 1

[David Goodmanson]

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))