Finite Difference Method & ODEs

FIZ353 - Numerical Analysis | 31/12/2020

Emre S. Tasci emre.tasci@hacettepe.edu.tr

At the end of our lecture on Interpolation, as a bonus, we talked about the finite difference method and employed it to solve an ordinary differential equation with respect to the given boundary conditions.

Taylor series expansion of the first derivative

Let's expand $f(x_i)$ to reach $f(x_{i+1})$ with $\Delta x \equiv x_{i+1} - x$:

$$f(x_{i+1}) = f(x_i)+\frac {f'(x_i)}{1!} (\Delta x)+ \frac{f''(x_i)}{2!} (\Delta x)^2+\frac{f'''(x_i)}{3!}(\Delta x)^3+ \cdots\quad(1)$$

re-arranging, we get:

$$f'(x_i) = \frac{f(x_{i+1})-f(x_{i})}{\Delta x} -\frac{f''(x_i)}{2!} \Delta x-\frac{f'''(x_i)}{3!}(\Delta x)^2+ \cdots$$

and if we expand $f(x_i)$ to reach $f(x_{i-1})$, this time, assuming we are maintaining the same distance length between consecutive 'steps', the difference $x_{i-1} - x_{i}$ will be $-\Delta x$

If this is not obvious, consider $x{i+1} = 6$, $xi = 5$ and $x{i-1}=4$ points: then $x_{i+1} - xi = 6-5=1=\Delta x$ and $x{i-1} - xi =4 - 5=-1=-\Delta x$.

so: $$f(x_{i-1}) = f(x_i)-\frac {f'(x_i)}{1!} (\Delta x)+ \frac{f''(x_i)}{2!} (\Delta x)^2-\frac{f'''(x_i)}{3!}(\Delta x)^3+ \cdots\quad(2)$$

Subtracting the second from the first and then isolating $f'(x_i)$ yields:

$$f'(x_i)= \frac{f(x_{i+1})-f(x_{i-1})}{2\Delta x}-\frac{f'''(x_{i})}{6}(\Delta x)^2+\cdots\quad(3)$$

The errors I've made in my life calculations...

Let's ponder on the order of the errors we suffer when we only take the first terms and ignore the rest, with respect to the three approaches above (namely, the "forward finite difference", the "backward finite difference" and the "centered finite difference" approaches -- for details, the reader is once again reffered to the bonus section of the interpolation lecture)

Take $f(x) = \sin(x)$ as our function. Let's say we want to calculate it's derivative at $x_i=\tfrac{\pi}{3}$ and $\Delta x = 0.5$

Forward Finite Difference Approximation From the (1) equation, taking upto the 1st term yields: $$f_f'(x_i) =\frac{f(x_{i+1})-f(x_{i})}{\Delta x}\underbrace{ -\frac{f''(x_i)}{2!} \Delta x-\frac{f'''(x_i)}{3!}(\Delta x)^2+ \cdots}_{\text{ignored part}} = \frac{f(x_{i+1}) - f(x_i)}{\Delta x} + O(\Delta x)$$

$$\Rightarrow\boxed{f_f'(x_i) \approx \frac{f(x_{i+1}) - f(x_i)}{\Delta x}}$$

Thus:

$$f_f'(\tfrac{\pi}{3}) \approx \frac{f(\tfrac{\pi}{3}+0.5)-f(0)}{0.5}$$

It doesn't look very promising. Let's see how the error value changes as we shorten our $\Delta x$ and while we're doing it, let's also calculate the first term approximations for the backward finite difference and the centered finite difference methods as well:

Backward Finite Difference Approximation From the (2) equation, taking upto the 1st term yields:

$$\boxed{f_b'(x_i) \approx \frac{f(x_{i}) - f(x_{i-1})}{\Delta x}}$$

Centered Finite Difference Approximation Directly from the (3) equation, taking upto the 1st term we have:

$$\boxed{f_c'(x_i)\approx \frac{f(x_{i+1})-f(x_{i-1})}{2\Delta x}}$$

It is obvious that the error of the centered finite difference approximation is much less compared to the forward & backward finite difference approximations. In order to see why is that, consider the higher terms (2nd and afterwards) that we are ignoring in our calculations: For the forward & backward methods, they start from the order of $\Delta x$, however for the centered method, the ignored part is in the order of $\Delta x$.

2nd order derivative approximation:

Once again, in our quick introduction to the topic last lecture, we had gone around and derived in a fast-forward way an identity/approximation for the 2nd order derivative using logic and the already found 1st order derivatives as:

$$\frac{\text{d}^2 f(x)}{\text{d}x^2}\biggr\vert_{x_i} \approx \frac{\frac{f(x_{i+1})-f(x_{i})}{h}-\frac{f(x_{i})-f(x_{i-1})}{h}}{h}=\frac{f(x_{i+1})-2f(x_{i})+f(x_{i-1})}{h^2}$$

Now, let's redrive it more gracefully. First remember our equations for the Taylor expansions of $f(x_{i+1})$ and $ f(x_{i-1})$:

$$f(x_{i+1}) = f(x_i)+\frac {f'(x_i)}{1!} (\Delta x)+ \frac{f''(x_i)}{2!} (\Delta x)^2+\frac{f'''(x_i)}{3!}(\Delta x)^3+ \cdots\quad(1)$$$$f(x_{i-1}) = f(x_i)-\frac {f'(x_i)}{1!} (\Delta x)+ \frac{f''(x_i)}{2!} (\Delta x)^2-\frac{f'''(x_i)}{3!}(\Delta x)^3+ \cdots\quad(2)$$

Then we ignore the 3rd and higher derivates and isolate the 2nd derivatives to the left hand side:

$$(1)\Rightarrow f''(x_i) \approx \frac{2!}{(\Delta x)^2} \left[f(x_{i+1}) - f(x_i) - \frac {f'(x_i)}{1!} (\Delta x)\right]\quad(1')$$$$(2)\Rightarrow f''(x_i) \approx \frac{2!}{(\Delta x)^2} \left[- f(x_i) + f(x_{i-1}) + \frac {f'(x_i)}{1!} (\Delta x)\right]\quad(2')$$

Adding them side by side will yield:

$$(1')+(2')\Rightarrow 2f''(x_i) \approx \frac{2!}{(\Delta x)^2} \left[f(x_{i+1}) - 2f(x_i) + f(x_{i-1})\right]$$$$\require{cancel}\cancel{2}f''(x_i) \approx \frac{\cancel{2!}}{(\Delta x)^2} \left[f(x_{i+1}) - 2f(x_i) + f(x_{i-1})\right]$$$$\boxed{f''(x_i) \approx \frac{f(x_{i+1}) - 2f(x_i) + f(x_{i-1})}{(\Delta x)^2}}$$

A Question and a Comment!

Question: Can you figure out the order of error of this above approximation? ;)

Comment: Boy, these derivations eerily look like the workings of the Newton polynomials (of the interpolation lecture), innit? ;;;)

Definition

An ordinary differential equation (ODE) is an equation that includes (usually linear) combinations of a function's derivatives with respect to its parameter. The "ordinary" indicates that, the function depends only on one parameter (as opposed to partial derivatives).

$$a_0(x)y +a_1(x)y' + a_2(x)y'' +\cdots +a_n(x)y^{(n)}+b(x)=0$$

Our very first ODE

Here's an equation you might be familiar with, i.e., the mass-spring system! From Hooke's Law, we have:

$$ F = -kx$$

and Hooke's nemesis, Sir Isaac Newton tells us that:

$$ F = m \ddot{x}$$

and therefore, combining these two by eliminating F, we get:

$$m\ddot x = -k x\rightarrow \ddot x(t) = -\frac{k}{m} x(t)$$

(I've explicitly indicated the time dependence of position)

and giving $\tfrac{k}{m}$ a fancy name like "$\omega^2$", we have:

$$\boxed {\ddot x(t) = -\omega^2 x(t)}\quad\left(\omega\equiv\sqrt{\frac{k}{m}}\right)$$

As opposed to the classic opening to the analitic solution that is "let's assume $x(t) =x_0 \cos{(\omega t + \varphi)}$, we will actually solve it without let's assuming anything:

image-2.png

We are given these data:

$$ k=1\text{ N/m}\\ m = 1\text{ kg}\\ x_0 = 10\text{ cm} = 0.1\text{ m}\\ v_0 = 0$$

So:

$$\omega = \sqrt{\frac{k}{m}} = 1\text{ s}^{-1}\\ x(t=0) = 0.1\text{ m}\\ \dot x(t=0) = 0$$

We write our approximation for the 2nd derivative:

$$f''(x_i) \approx\frac{f(x_{i+1}) - 2f(x_i) + f(x_{i-1})}{(\Delta x)^2}$$

In our problem, this translates to (I'm now switching to equality sign even though it is approximation for practicality): $$\ddot x_0 = \frac{x_1 - 2x_0 + x_{-1}}{(\Delta t)^2}$$

The $0$ indice denotes the value at $t=0$, $1$ indice the value at $t=\Delta t$ and $-1$ denotes the value at $t=-\Delta t$.

Plugging this into the $\ddot x(t) = -\omega^2 x(t)$ equation yields:

$$\frac{x_1 - 2x_0 + x_{-1}}{(\Delta t)^2} = -\omega^2 x_0$$

Looking at this equation, we realize that, if we knew the positions $x_0$ and $x_{-1}$ we could recover $x_1$ and increase the indice & re-iterate it as:

$$\frac{x_2 - 2x_1 + x_{0}}{(\Delta t)^2} = -\omega^2 x_1$$

and it would continue on and on...

Let's focus on our initial equation once again:

$$\frac{x_1 - 2x_0 + x_{-1}}{(\Delta t)^2} = -\omega^2 x_0$$

We know $x_0$ but not $x_{-1}$. Can't we use our heads?...

Let's put our physics knowledge in front of math: We have pulled the spring up to $x_0 = 0.1\text{ m}$ and then let it go. From physics we know that, this is the maximum displacement our spring-mass will achieve. We can also deduce the fact that, the movement will be symmetrical with respect to the maximum point, e.g., the particle will be at the same position for $t=-\Delta t$ and $t=\Delta t$, meaning that $x_1 = x_{-1}$. Substituting $x_1$ as $x_{-1}$, our equation becomes:

$$\begin{align*}\frac{x_1 - 2x_0 + x_{-1}}{(\Delta t)^2} &= -\omega^2 x_0\\ \frac{x_1 - 2x_0 + x_{1}}{(\Delta t)^2} &= -\omega^2 x_0\\ \frac{2x_1 - 2x_0}{(\Delta t)^2} &= -\omega^2 x_0 \end{align*}$$

Re-arranging yields:

$$x_1 = 2x_0 - \omega^2 {\Delta t}^2 x_0 = \left(1-\frac{\omega^2 {\Delta t}^2}{2}\right)x_0$$

But don't forget that, this last equation we derived is only correct for $x_{-1}$ and $x_1$ -- it is not a general relation as:

$$\frac{x_{i+1} - 2x_i + x_{i-1}}{(\Delta t)^2} = -\omega^2 x_i$$

which corresponds to:

$$ x_{i+1} = \underbrace{\left(2 -\omega^2{\Delta t}^2\right)}_{\zeta} x_i - x_{i-1}$$

We identified $\left(2 -\omega^2{\Delta t}^2\right)$ as $\zeta$ because it doesn't change its value and it would be a bad example to recalculate it at every iteration, so:

$$ x_{i+1} =\zeta x_i - x_{i-1}$$

Now that we have the get going, let's give it a try!

Well, we got something to ourselves, so let's see if this is indeed what we wanted by having plotting it:

We can see the period being something around 6 seconds, let's also check it out by first creating a filter that will be true if the position of that row is equal to the position of the first row (i.e., 0.1) and then, apply this to our calculated data:

Only the first entry (being equal to itself -- duh!)??? Why can it be? Oh, right, that's because these are NUMERICAL VALUES OBTAINED FROM NUMERICAL VALUES! So there's a very low probability that we will have a definite "0.1" but they will rather be something like "0.100000000000001" or "0.099999999999999998" - always keep this in mind when filtering!

So instead of seeking those entries with positions exactly equal to that of the first row's position, we go for very close values -- it would be a miracle if our designated $\Delta t$ would exactly step on the period of the oscillation! Therefore we will take this as an "equality":

$$|x_i - x_0| < 10^{-6} \rightarrow x_i \approx x_0$$

and here is how we write it in code:

The period seems to be 6.27 seconds. Let's calculate it using the relation between the angular frequency $\omega$ and the period:

$$ T = \frac{2\pi}{\omega}$$

A Question and a Hint

Question: We could actually do much better than taking the first period we see... what was it that we could do?...

Hint:

When solving this case, we actually kind of 'cheated' somewhere at the start (do you see where it was?)...

Yes, it was the physics thing ($x1 = x{-1}) - we used our physics head to get over the difficult part. Could we have done it without resorting to physical interpretation, just treating it as a math problem? Yes we could. First the implementation, then the explanation:

Actually, for the sake of the discussion, let's focus on the 1st period:

Now we can talk about what we did: in the second approach, we just assigned the same value of $x_0$ to $x_1$ and moved on from there. Does this make sense? As long as our incrementation $\Delta t$ is small enough.

Let's follow the same approach but this time taking $\Delta t = 0.1$:

Now we start to spot the difference. To bluntly put it:

In summary, when we take $\Delta t$ sufficiently small, under a reasonable force(/acceleration) the system will not differ much from its state in the previous iteration. This is exactly what we are exploiting by using small differences. Its downside is the computational cost: the smaller the differences are, the greater the number of calculation points and hence the longer the calculation time. There are of course a couple of optimization techniques such as updating the incrementation size with respect to the slope of the function at the current point, etc.

Systematical Approaches

Now that we have wildly entered into the scene, it's time for more systematical approaches

Euler's method

This method is actually the proper treatment to what we have done so far, only with a limitation to 1st order derivatives(*) and with the initial value given: We approximate the derivatives with slopes, such that:

$$\text{Given: }y'(x) = f(x,y(x)),\;y(x_0) = y_0\\ \Rightarrow y_{n+1} = y_n + \Delta x f(x_n,y_n)$$

where $f(x_n,y_n)$ is called as the incrementation function $\varphi$.

(*) Even though Euler's method works for 1st order derivatives, any ODE of order N can be re-written as a system of 1st order ODE equations.

Example

Use Euler's method to solve

$$y' = 4e^{0.8x} - 0.5y$$

from $x=0$ to $x=4$ with a step size of $\Delta x = 1$

with the given initial condition: $y(x=0)=2$.

Compare your results with the analytical solution:

$$y(x) = e^{-0.5x}\left(3.07692e^{1.3x}-1.07692\right)$$

image.png

Solution $$y(1) = y(0) + 1.f(x_0=0,y_0=2)\\ f(0,2) = 4e^0-0.5\times2=3\\ \rightarrow y(1) = 2 +1\times3=5$$

$$y(2) = y(1) + 1.f(x_1=1,y_1=5)\\ f(1,5) = 4e^{0.8\times 1}-0.5\times5=6.40216\\ \rightarrow y(2) = 5 +1\times6.40216=11.40216\\ \vdots$$

Let's decrease the step size:

Runge - Kutta Methods

Runge-Kutta (RK) methods can be thought as a higher order version of the Euler's method. The incrementation factor $\varphi$ is defined as:

$$\varphi = a_1 k_1+a_2 k_2+\dots+a_n k_n$$

where $a$'s are constants and the $k$'s are defined as:

$$k_1 = f(x_i,y_i)\\ k_2 = f(x_i+p_1\Delta x,y_i+q_{11}k_1\Delta x)\\ k_3 = f(x_i+p_2\Delta x,y_i+q_{21}k_1\Delta x + q_{22}k_2\Delta x)\\ \vdots\\ k_n = f(x_i+p_{n-1}\Delta x,y_i + q_{n-1,1}k_1 \Delta x+ q_{n-1,2}k_2 \Delta x+\dots+q_{n-1,n-1}k_{n-1} \Delta x)$$

where the $p$'s and $q$'s are constants.

(for $n=1$ RK converges to Euler's method)

Classic Runge-Kutta Method

Given $\frac{dy}{dx} = f(x, y), \quad y(x_0) = y_0$

The most employed RK method is the 4th order RK approach with the following fix:

$$\begin{align} x_{n+1} &= x_n + \Delta x \\ y_{n+1} &= y_n + \frac{1}{6}\Delta x\left(k_1 + 2k_2 + 2k_3 + k_4 \right) \end{align}$$

with:

$$\begin{align} k_1 &= \ f(x_n, y_n), \\ k_2 &= \ f\left(x_n + \frac{\Delta x}{2}, y_n + \Delta x\frac{k_1}{2}\right), \\ k_3 &= \ f\left(x_n + \frac{\Delta x}{2}, y_n + \Delta x\frac{k_2}{2}\right), \\ k_4 &= \ f\left(x_n + \Delta x, y_n + \Delta xk_3\right). \end{align}$$

image.png
[Source: Chapra]

image-2.png
Source: Wikipedia, By HilberTraum - Own work, CC BY-SA 4.0

Example

Use RK method to solve

$$y' = 4e^{0.8x} - 0.5y$$

from $x=0$ to $x=4$ with a step size of $\Delta x = 1$

with the given initial condition: $y(x=0)=2$.

Compare your results with the analytical solution:

$$y(x) = e^{-0.5x}\left(3.07692e^{1.3x}-1.07692\right)$$

Example: Bungee Jumper

A bungee jumper's acceleration, with the drag force (which is proportional to $v^2$) taken into account is given by:

$$a = g - \frac{c_d}{m}v^2$$

where $g=9.81$ m/s2, $m=68.1$ kg and the drag coefficient is $c_d = 0.25$ kg/m.

Solve for the position & velocity with respect to time, assuming that $x(t=0)=v(t=0) = 0$ and take $\Delta t =h= 2$ s increments.

Compare your findings with the analytical solutions: $$x(t)=\frac{m}{c_d}\ln{\left[\cosh\left({\sqrt{\frac{g c_d}{m}}t}\right)\right]}$$

$$v(t)=\sqrt{\frac{gm}{c_d}}\tanh\left({\sqrt{\frac{g c_d}{m}}t}\right)$$

Solution $$\frac{\text{d}x}{\text{d}t}=v$$ $$\frac{\text{d}v}{\text{d}t}=a=g - \frac{c_d}{m}v^2\\ \downarrow$$

$$\frac{\text{d}x}{\text{d}t}=f_1(t,x,v)=v\rightarrow k_{1,1}=f_1(0,0,0)=0\\ \frac{\text{d}v}{\text{d}t}=f_2(t,x,v)=g- \frac{c_d}{m}v^2\rightarrow k_{1,2}=f_2(0,0,0)=9.81 - \frac{0.25}{68.1}(0)^2 = 9.81$$$$x(1) = x(0) + k_{1,1}\frac{h}{2}=0+0\frac{2}{2} = 0\\ v(1) = v(0)++ k_{1,2}\frac{h}{2}=0+9.81\frac{2}{2} = 9.81$$

Now that we have the $k_1$'s, we can proceed to calculate the $k_2$'s:

$$k_2 = f(t_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1 h)\\ k_{2,1}=f_1(1,0,9.81)=9.8100\\ k_{2,2}=f_2(1,0,9.81)=9.4567\\ x(1) = x(0) + k_{2,1}\frac{h}{2}=0+9.8100\frac{2}{2}=9.8100\\ v(1) = v(0) + k_{2,2}\frac{h}{2}=0+9.4567\frac{2}{2}=9.4567 $$

Then comes the $k_3$'s:

$$k_3 = f(t_i+\frac{1}{2}h,y_i+\frac{1}{2}k_2 h)\\ k_{3,1}=f_1(1,9.8100,9.4567)=9.4567\\ k_{3,2}=f_2(1,9.8100,9.4567)=9.4817\\ x(2) = x(0) + k_{3,1}h=0+9.4567\times2=18.9134\\ v(2) = v(0) + k_{3,2}h=0+9.4817\times2=18.9634 $$

Finally, the $k_4$'s:

$$k_4=f(t_i+h,y_i+k_3 h)\\ k_{4,1}=f_1(2,18.9134,18.9634)=18.9634\\ k_{4,2}=f_2(2,18.9134,18.9634)=8.4808$$

combining them all, we have:

$$\boxed{x(2)=0+\frac{1}{6}\left[0+2\left(9.8100+9.4567\right)+18.9634\right]\times 2=19.1656\\ v(2)=0+\frac{1}{6}\left[9.8100+2\left(9.4567+9.4817\right)+8.4898\right]\times 2=18.7256}$$

image.png
[Source: Chapra]

RK45 function

In 1980, combining the fourth and fifth order RK formulas, Dormand and Prince came up with a very efficient and widely used RK approach, implemented to libraries as RK45.

In Python, this function is implemented as scipy.integrate.RK45

solve_ivp function

Yet another function that acts as a compilation of solvers (with RK45 being the default) is the scipy.integrate.solve_ivp function.

This function is much easier to implement and use as evaluation points can be specified via the _tspan parameter and unlike the above RK45 function, it isn't a solver object, therefore there isn't any requirement to advance via the step() command.

Setting the dense_output option to True, we get the .sol method that acts as a function and evaluates the solution at the given t.

now that we have the $v$ values, we can also calculate the $x$ values as well. The ODE we need to solve is straightforward:

$$\dot x = v$$

and to calculate the $v$ for the control points, we are implement a function that just feeds these to the solver! 8)

As a final exercise, let's also have the real values calculated and re-drive the results table that was quoted above:

$$x(t)=\frac{m}{c_d}\ln{\left[\cosh\left({\sqrt{\frac{g c_d}{m}}t}\right)\right]}$$$$v(t)=\sqrt{\frac{gm}{c_d}}\tanh\left({\sqrt{\frac{g c_d}{m}}t}\right)$$

References