Python code for Second Order ODE Initial Value problem using Finite Difference methods

64 Views Asked by At

I am trying to solve the following ODE $$\begin{aligned} & y^{\prime \prime}+2 y^{\prime}+y=0 \\ & y(0)=2 \quad y^{\prime}(0)=-1 \\ & 0 \leqslant x \leqslant 1 \end{aligned}$$ Using central FD, I can write $$\frac{y_{i+1}-2 y_i+y_{i-1}}{h^2}+2 \frac{y_{i+1}-y_{i-1}}{2 h}+y_i=0$$ which I have simplied to $$(1+h) y_{i+1}+\left(h^2-2\right) y_i+(1-h) y_{i-1}=0$$ I am struggling to fit the second boundary condition into the set of resulting equations. Using central FD, I think I can derive $$\begin{aligned} & \frac{y_2-y_0}{2 h}=-1 \\ & y_0=y_2+2 h\end{aligned}$$ However, I am not sure how to make use of this condition. Can you please guide me ?

The analytic solution to the ODE is given as $y(x)=2 e^{-x}+x e^{-x}$

Thanks

Based on whpowell96's suggestion, I will create the system of equations as follows. I have removed the python code for now. Does this look correct? Thank you.

Assuming N = 10

$$ \begin{aligned} & h=0.125 , h^2=0.016, 1-h=0.875 \\ & 1.125 y_2-1.984 y_1=-1.75 \\ & 1.125 y_3-1.984 y_2+0.875 y_1=0 \\ & 1.125 y_{11}-1.984 y_{10}+0.875 y_9=0 \end{aligned} $$

AND from $B C$ $$ 1 y_2+0 y_1-1 y_0=-2 h $$

$$\begin{bmatrix} -1.984 & 1.125 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0.875 & -1.984 & 1.125 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0.875 & -1.984 & 1.125 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0.875 & -1.984 & 1.125 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0.875 & -1.984 & 1.125 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0.875 & -1.984 & 1.125 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0.875 & -1.984 & 1.125 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0.875 & -1.984 & 1.125 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.875 & -1.984 & 1.125 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.875 & -1.984 & 1.125\\ -1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \\ y_6 \\ y_7 \\ y_8 \\ y_9 \\ y_{10} \\ y_{11} \\ \end{bmatrix} = \begin{bmatrix} -1.75\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ -0.25\\ \end{bmatrix}$$

1

There are 1 best solutions below

0
On

The way to do this that achieves the best accuracy is via the ghost point method. Label the nodes such that $y_0$ is at $x=0$. The discretized boundary conditions read $$ \begin{aligned} y_0 &= 2 \\ \frac{y_{1} - y_{-1}}{2h} &= -1. \end{aligned} $$ Here we have introduced an extra variable, $y_{-1}$ which corresponds to the function value one node to the left of the domain. Since the function is not defined there, it is called a ghost node. We won't actually evaluate it, but will merely use it to inform our algebra.

Now consider the difference equation at $x=0$: $$ (1+h)y_1 + (h^2-2)y_0 + (1-h)y_{-1} = 0. $$ Again, we are using the ghost point $y_{-1}$. However, notice that our boundary condition can be rearranged to $y_{-1} = 2h + y_1$. Substituting this and $y_0=2$ into the difference equation yields $$ (1+h)y_1 + 2(h^2-2) + (1-h)(2h+y_1)=0, $$ which is a linear equation for $y_1$. Since the boundary conditions was discretized using a second order scheme, we actually retain second order accuracy of the finite difference scheme, which is lost if you approximate the boundary condition with a simple one-sided difference.