Trying to model a simple second order ODE

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.

• The general approach to solving second- or higher-order ODEs is to first convert them into a system of first-order equations. Dec 11 '20 at 2:24
• this example may or may not be helpful in some way
– uhoh
Dec 11 '20 at 2:51
• Initial conditions are not clear to me: the first one is $y(0)=0$, but it seems the other one is $y(h)=-h^2/2 - h$, instead of $y'(0)$.
– VoB
Dec 11 '20 at 11:00
• @VoB : The method that is implemented is a linear multi-step methods with 2 memory positions. Which means that the first step will compute $u_2$ from $u_0$ and $u_1$. $u_1$ has to be determined by some other means from $u(0)$ and $u'(0)$, the error order of $u_1$ has to be at least as high as the method order. It is not unreasonable to take the exact value there for test purposes. Dec 11 '20 at 14:26
• @LutzLehmann Yes, thanks. I wrote Backward Euler in my answer just because OP mentioned BE in the question, but I agree with you
– VoB
Dec 11 '20 at 14:54

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$$ Due to the first order approximation of the first derivative, the coefficients $$C$$, $$D$$ will grow like $$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.

• Thank you so much. This is due to fact that I mix the first and second orders derivatives or just to the poor approximation? Many thanks!! Dec 11 '20 at 13:10
• Not really. If you change it to the central difference quotient, this increases the truncation error to second order. But that only means that the coefficients $C,D$ move from zero like $h^2x$. This removes about 2 orders of magnitude from the scale, so from 41 to 39 or 38 in the power of 10. It might happen that the direction changes, so that the plot goes to the negative large numbers. Dec 11 '20 at 13:27
• Thank you @Lutz Lehmann Dec 11 '20 at 13:33
• There is no way to avoid that error. The $y$-Lipschitz constant is trivially $1$, the floating-point errors thus get magnified with factor $e^{Lx}=e^x$ and the factor to translate the local error into the global error is generally $\frac{e^{Lx}-1}{Lh}$. So independent of the method of error $p$ you will get a global error of $O(e^{x}(h^p+\mu))$. So to integrate over longer intervals than $T=35$ you need a multi-precision data type with a correspondingly smaller "machine constant" $\mu$, and need to adapt stepsize and order to fit. Dec 11 '20 at 17:38
• In practical situations one of course avoids constructions that have exponential or faster explosions. This typically results in bounded solutions with transients towards a low-dimensional attractor manifold. In these circumstances the accumulated errors tend to have an oscillating sign and thus cancel. It sometimes seems that the largest error source is a time dilation effect, which will of course lead to large errors at larger times without changing much of the qualitative behavior. Dec 11 '20 at 17:43

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

$$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.

• Thank you so much. Could you maybe cite a reference in wich there are explanations why could not I use higher orders discretizations along with others orders? Is there something wrong with scale perhaps? Many thanks!! Dec 11 '20 at 13:09
• Any good numerical analysis book will discuss this point: keep in mind the mantra "higher order requires higher regularity". Think about finite differences: you can obtain pretty good approximations (i.e. increasing the order) as long as your function is regular enough: that's because you are relying on Taylor expansions @Still_waters
– VoB
Dec 11 '20 at 13:52
• It appears that this solver follows the polynomial path as long as the floating point errors (of size $\mu$) are not amplified into visibility, that is, until $\mu e^x\sim 1$, which gives $x\approx 37$. After that the exponential regime takes over and quickly diverges from the exact solution. Dec 11 '20 at 15:40
• @LutzLehmann Sure, that's the problem with long time integration for simple methods... but I thought the OP's question was about how to solve it. Btw, thanks for pointing that out! :-)
– VoB
Dec 11 '20 at 17:53
• @Still_waters just an English note: I had no idea what 'principiant' meant until I googled it, I think 'I'm a novice' will better communicate what you're trying to say Dec 11 '20 at 18:11