5
$\begingroup$

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

enter image description here

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.

$\endgroup$
12
  • 4
    $\begingroup$ The general approach to solving second- or higher-order ODEs is to first convert them into a system of first-order equations. $\endgroup$ Dec 11, 2020 at 2:24
  • 1
    $\begingroup$ this example may or may not be helpful in some way $\endgroup$
    – uhoh
    Dec 11, 2020 at 2:51
  • 1
    $\begingroup$ 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)$. $\endgroup$
    – VoB
    Dec 11, 2020 at 11:00
  • 2
    $\begingroup$ @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. $\endgroup$ Dec 11, 2020 at 14:26
  • 2
    $\begingroup$ @Still_waters There is nothing wrong *in principle" with using higher order discretizations. What you did wrong was to pick some discretizations apparently at random, without doing any form of error analysis on your method. That is almost never going to work. If you want to do some more "random" experiments, try using forwards, central, and backward approximations for both derivatives (which will give you 9 different methods to try). $\endgroup$
    – alephzero
    Dec 11, 2020 at 22:59

2 Answers 2

5
$\begingroup$

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.

relative error, w/o modifications

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

relative error, Kahan

$\endgroup$
17
  • 2
    $\begingroup$ 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. $\endgroup$ Dec 11, 2020 at 13:27
  • 3
    $\begingroup$ 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. $\endgroup$ Dec 11, 2020 at 17:38
  • 2
    $\begingroup$ 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. $\endgroup$ Dec 11, 2020 at 17:43
  • 1
    $\begingroup$ @LutzLehmann Why is the $y$ Lipschitz constant 1? Did you reduce to first order system, so the r.h.s. is $f(y,y',x)=[y',y'+x]$ and hence the Jacobian is $J_f = $\begin{bmatrix} 0 & 1 \\ 0 & 1 \end{bmatrix} which has inf norm equal to 1? $\endgroup$
    – Vefhug
    Dec 30, 2020 at 20:17
  • 1
    $\begingroup$ @Vefhug : You are right, "trivially $1$" is only true if the order is reduced via $v=y'$. For the second order system one needs an adapted norm, as you computed the max norm giving the row-sum norm is suitable here. In other norms one could get other constants, but I think all give $1\le L\le 2$, so the principal argument about the error growth remains valid. $\endgroup$ Dec 30, 2020 at 20:54
7
$\begingroup$

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

enter image description here

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.

$\endgroup$
7
  • $\begingroup$ 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!! $\endgroup$ Dec 11, 2020 at 13:09
  • 4
    $\begingroup$ 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 $\endgroup$
    – VoB
    Dec 11, 2020 at 13:52
  • 2
    $\begingroup$ 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. $\endgroup$ Dec 11, 2020 at 15:40
  • 2
    $\begingroup$ @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! :-) $\endgroup$
    – VoB
    Dec 11, 2020 at 17:53
  • 2
    $\begingroup$ @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 $\endgroup$
    – llama
    Dec 11, 2020 at 18:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.