Trying to model a simple second order ODE

Question

I am studying some computational methods and I am trying to program simples equations to understand how the methods work... Particularly, I am trying to understand how multiorders ODE's work.

I've tried to program the ODE $y''=y'+x$, with conditions $y(0)=0$ and $y(h)=-h^2/2-h$, ie, the function $y=-x^2/2-x$.

So, I've used central diferences on second order and backward Euler on first order, getting

$$\dfrac{u_{i+1}-2u_i+u_{i-1}}{h^2}=\dfrac{u_i-u_{i-1}}{h}+x_{i+1}$$

$$u_{i+1}=u_i(h+2)+u_{i-1}(-h-1)+h^2x_{i+1}$$.

However, the result is terrible.

h=0.05;
max=100;
np=max/h;
u=[];
y=[];
u(1)=0;
u(2)=-h^2/2-h;
y(1)=0;
y(2)=-h^2/2-h;
x=[0:h:max];
for i=3:np+1;
    u(i)=u(i-1)*(h+2)+u(i-2)*(-h-1)+h^2*x(i);
    y(i)=-(x(i))^2/2-(x(i));
end
close all
plot(x,u)
hold on
plot(x,y,'r')

I know there is a lot of better methods etc, but I'd like to know why this terrible result. I mean, I didn't expect this impressive error. I am thankful in advance.

Answer 1

The result you see is due to the long integration interval and the accumulation of truncation errors. Each step with its truncation error will switch to a slightly different exact solution. Now the exact solutions are in general $$ y(x)=-\frac12x^2-x+C+De^x, ~~~y(0)=0\implies C+D=0. $$ The numerical solution will move trough a sequence of exact solutions that are determined by the value pairs $u_n,u_{n+1}$. Due to the first order approximation of the first derivative, the coefficients $C$, $D$ will grow like $h^2n=hx$. This gives a value of $5$ at $x=100$ for $h=0.05$. But the biggest contribution is due to the exponential, as $e^{100}=2.688117·10^{43}$. This is indeed the magnitude you see in the plot.

To get a more appreciable measure, the relative error is smaller $0.1$ up to $x=4.1$, if the iteration formula is corrected to use x(i-1) to compute u(i).

---

On a second look

Already setting the value $u_2$ to the exact value $y(h)$ introduces a non-zero exponential factor into the numerical solution $$ u_n=Ax_n^2+Bx_n+C((1+h)^{n-1}-1),~~n=1,2,... $$ as the polynomial solutions of differential and difference equation are different. The values $A,B$ for the particular polynomial solution are determined by insertion into the difference equation, \begin{align} \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}&=\frac{u_i-u_{i-1}}{h}+x_i\\ 2A&=(A(2x_i-h)+B)+x_i\\ 2A &=-Ah +B \\ 0&=2A+1\\ \hline A=-1/2,&~~B=-1-h/2 \end{align}

To get $C=0$ in the exact solution of the difference equation, the second starter value has to be chosen as $u_2=Ah^2+Bh=-h^2-h$.

With the values from the ODE solution $u_1=0$, $u_2=-h-h^2/2$ we get for the general solution \begin{align} u_n&=-x_n^2/2-(1+h/2)x_n+C((1+h)^{n-1}-1) \\ -h^2/2-h=u_2&=-h^2-h+Ch \\[.5em] h/2&=C \end{align} Thus $C=h/2$ for a total error of $$ u_{n+1}-y(x_{n+1})=h/2·((1+h)^n-1-nh). $$ Also in this theoretical estimate the relative error is smaller than $0.1$ for $x\le 4.1$, with the absolute error smaller $0.1$ up to $x=2$.

With the corrected starter values the plot for the actual numerical computation gives a relative error of $0.1$ or smaller up to $x=34.3$. In the error plot the two phases where the error behaves like the theoretical $h/2/(1+x)$ and then switches to exponential growth are clearly recognizable.

Using Kahan summation for the accumulation of the step updates for an effective computation in quadruple precision gives a relative error smaller $0.1$ up to $x=44.0$.

Answer 2

I got your point: you have the ODE $y''=y'+x$, and we can see that $y(x)=-\frac{x^2}{2} - x$ solves the problem. As you want to integrate it numerically, you want to set two initial conditions, so you imposed $y(0)=0$ (and that's fine) and also another one $y(h)$, which is wrong because it's not the required information you need to solve an ODE, mainly because it's not initial condition. What you need is $y'(0)=-1$

Given that, you can recast everything into a first order system (see @WolfgangBangerth's comment)

\begin{cases} y' = u \\ u' = u + x \\ y(0)=0 \\ u(0)=-1 \end{cases}

Now you can choose the time integration scheme you prefer. Since you mentioned Backward Euler, just set the vector valued function $f(x,Y) = [u,u+x]$, where $Y=[y,u]$ is your solution vector, and then solve $$Y_{n+1} = Y_n + h f(x_{n+1},Y_{n+1})$$ at each time step. Notice, however, that you can avoid to use a non-linear solver, as the r.h.s. may be written as $$A Y +b$$ where

$A= \begin{bmatrix} 0 & 1 \\ 0 & 1 \end{bmatrix}$ and $b= \begin{bmatrix} 0 \\ x \end{bmatrix}$

and hence your ODE becomes

$$Y' = A Y + b$$

and Backward Euler reads $$Y_{n+1} = Y_n + h(A Y_{n+1} + b(x_{n+1}))$$ i.e.

$$Y_{n+1} =(I_2 - hA)^{-1}(Y_n + h b(x_{n+1}))$$

This can be written in a simple Python snippet:

import numpy as np
import matplotlib.pyplot as plt

def sol(x):
    return -.5*x**2 - x

def b(x):
    return np.array([0,x])
    
A = np.array([[0,1],[0,1]])
Y0 = np.array([0.0,-1.0])
h = 0.01
T = 1.0
n = int(T/h)
x = np.linspace(0,T,n+1)
Y = np.zeros([2,n+1])
Y[:,0] = Y0.copy()
I = np.eye(2)
error = np.zeros(n)
for i in range (0,n):
    Y[:,i+1] = np.linalg.solve(I-h*A, Y[:,i]+ h*b(x[i+1]))

plt.plot(x,Y[0,:],'o',markerfacecolor='None')
plt.plot(x,sol(x),'-r')
plt.show()

which reproduce the solution correctly

Notice that I wrote a snippet for this particular case, with this particular r.h.s. just because it's linear! Otherwise, a Newton method would have been employed to solve the non-linear system at each time step. This latter approach can be applied wherever your r.h.s. $f$ is.