Numerical Solutions of Partial Differential Equations

First we solve an ordinary differential equaiton

$$\frac{d^2 u}{dx^2}-12 xu =-1$$ with $0 < x < 1$, along with the boundary conditions $u(0)=1$ and $u(1)=-1.$.

To solve this problem numerically, we divide the interval $[0,1]$ into n equal subintervals with length $\Delta x =\frac{x_n-x_0}{n}$. For our case $\Delta x = \frac{1-0}{n}=\frac{1}{n}$. The set of partitioning points or mesh points are given as

$$x_0=0$$ $$x_1=x_0+\Delta x=0+\frac{1}{n}=\frac{1}{n}$$ $$x_2 = x_0+2 \Delta x = 2 \Delta x =\frac{2}{n}$$ $$\vdots$$ $$x_i= x_0+i \Delta x =0+\frac{i}{n}$$ $$x_n=1$$ We approximate the terms of the differential equation as follows: $$u(x)\rightarrow u_i$$ $$\frac{du}{dx}\rightarrow \frac{u_{i+1}-u_{i-1}}{2 \Delta x}$$ $$\frac{d^2u}{dx^2}\rightarrow \frac{u_{i+1}-2u_i+u_{i-1}}{(\Delta x)^2}$$ $$u(x_0)\rightarrow u_0$$ $$u(x_n)\rightarrow u_i$$ $$u_i\rightarrow u(x_i), i=0, 1, 2, \cdots, n$$ Then our original differential boundary value problem becomes the following n system of equations $$\frac{u_{i+1}-2u_i+u_{i-1}}{(\Delta x)^2} -12 x_i u_i =-1, \quad i=1, 2, \cdots, n-1$$ $$u_0=1, u_n=-1$$ For our purpose we use $n =5$. Then the above system becomes $$25(u_2-2u_1+u_0)-\frac{12}{5}u_1=-1$$ $$25(u_3-2u_2+u_1)-\frac{25}{4}u_1=-1$$ $$25(u_4-2u_3+u_2)-\frac{36}{5}u_1=-1$$ $$25(u_5-2u_4+u_3)-\frac{48}{5}u_1=-1$$ Using the boundary conditions $u_0=1, u_n=-1$ and simplifying, we get $$-52.4 u_1+25 u_2 = -26$$ $$25 u_1-54.8u_2+25u_3=-1$$ $$25 u_2 -57.2 u_3+25 u_4 =-1$$ $$25 u_3-59.6 u_4 =24$$ We can represent the above system of equations as a matrix equation $$A u = f$$ where, $$A=\begin{bmatrix} -52.4 & 25 & 0 & 0\\ 25 & -54.8 & 25 & 0\\ 0 & 25 & -57.2 & 24\\ 0 & 0 & 25 & -59.6 \end{bmatrix}$$ $$u =\begin{bmatrix} u_1 \\ u_2\\ u_3\\ u_4 \end{bmatrix}$$ and $$f= \begin{bmatrix} 26\\-1\\-1\\24\end{bmatrix}$$ The solution is $$u = A^{-1} f$$ The following Python program solves this equaiton.

---

Heat Equaiton

Consider the one dimensional heat equaiton

$$u_{t}=u_{xx}, \quad 0 < x < 1, 0 < t,$$ $$u(0,t)=0, u(1,t)=0, 0 < t,$$ $$u(x,0) =f(x), 0 < x < 1.$$

We use the mesh points

$${x_0, x_1, x_2, \cdots,x_i, \cdots, x_n}$$ along the spatial direction and the mesh points $${t_0, t_1, t_2, \cdots, t_m, \cdots}$$ along the time direction. We approximate the solution $u(x,t)$ by $u(x_i, t_m) =u_i(m)$. The partial derivative terms are approximated as $$\frac{\partial u}{\partial x}(x_i,t_m) \rightarrow \frac{u_{i+1}(m)-u_{i-1}(m)}{2 \Delta x}$$ $$\frac{\partial^2 u}{\partial x^2}(x_i,t_m) \rightarrow \frac{u_{i+1}(m)-2u_i(m)+u_{i-1}(m)}{(\Delta x)^2}$$ $$\frac{\partial u}{\partial t}(x_i,t_m) \rightarrow \frac{u_{i}(m+1)-u_{i}(m)}{\Delta t}$$ Then the above heat equaiton has the following discritized form: $$\frac{u_{i}(m+1)-u_{i}(m)}{\Delta t}= \frac{u_{i+1}(m)-2u_i(m)+u_{i-1}(m)}{(\Delta x)^2}, \quad i=1,2,\cdots,n-1, \quad m=0,1,2,\cdots$$ These equations are solved for $u_i(m+1)$ $$u_i(m+1) =r u_{i-1}(m)+(1-2r) u_i(m)+ru_{i+1}(m),$$ where $r =\Delta t / (\Delta x)^2$. This equaiton is valid only for $0 < r \leq 1/2$. This equation can be used to compute $u_i(m+1)$ from $u_i(m)$ at the preceding time level. As for example, the values of the $u_i$ at time level $m=1$ can be calculated by setting $m=0$ in the above equation as $$u_i(1)=ru_{i-1}(0)+(1-2r)u_i(0)+ru_{i+1}(0) \quad i=1,2,\cdots, n-1.$$ This process can be repeated to compute $u_i(m)$ for $m=2, 3, \cdots$.

For a special case, we take $\Delta x =1/4$ and $r=1/2$. This will give $\Delta t =1/32$. The equations at $i = 1, 2, 3$ become

$$u_1(m+1)=\frac{1}{2}(u_0(m)+u_2(m)),$$ $$u_2(m+1)=\frac{1}{2}(u_1(m)+u_3(m)),$$ $$u_3(m+1)=\frac{1}{2}(u_2(m)+u_4(m)).$$

The following Python program can be used to compute the numerical solutions of the heat equation.

---

Wave Equaiton

Consider the one dimensional wave equaiton

$$u_{tt}=u_{xx}, \quad 0 < x < 1, 0 < t,$$ $$u(0,t)=0, u(1,t)=0, 0 < t,$$ $$u(x,0) =f(x), 0 < x < 1.$$ $$\frac{\partial u}{\partial t}(x,0) =g(x), 0 < x < 1.$$

We use the mesh points

$${x_0, x_1, x_2, \cdots,x_i, \cdots, x_n}$$ along the spatial direction and the mesh points $${t_0, t_1, t_2, \cdots, t_m, \cdots}$$ along the time direction. We approximate the solution $u(x,t)$ by $u(x_i, t_m) =u_i(m)$. The partial derivative terms are approximated as $$\frac{\partial u}{\partial x}(x_i,t_m) \rightarrow \frac{u_{i+1}(m)-u_{i-1}(m)}{2 \Delta x}$$ $$\frac{\partial^2 u}{\partial x^2}(x_i,t_m) \rightarrow \frac{u_{i+1}(m)-2u_i(m)+u_{i-1}(m)}{(\Delta x)^2}$$ $$\frac{\partial u}{\partial t}(x_i,t_m) \rightarrow \frac{u_{i}(m+1)-u_{i}(m)}{\Delta t}$$ $$\frac{\partial^2 u}{\partial t^2}(x_i,t_m) \rightarrow \frac{u_i(m+1)-2u_i(m)+u_i(m-1)}{(\Delta t)^2}$$ Then the above wave equaiton has the following discritized form: $$\frac{u_i(m+1)-2u_i(m)+u_i(m-1)}{(\Delta t)^2}= \frac{u_{i+1}(m)-2u_i(m)+u_{i-1}(m)}{(\Delta x)^2}, \quad i=1,2,\cdots,n-1, \quad m=0,1,2,\cdots$$ After solving these equations for $u_i(m+1)$, we get $$u_i(m+1) =\rho^2 u_{i-1}(m)+2(1-\rho^2)u_i(m)+\rho^2 u_{i+1}(m)-u_i(m-1),$$ where $\rho =\Delta t / \Delta x$. For each time level, we need the informaiton for two previous levels. From the initial condition, we have $u_i(0)=f(x_i)$. To get the values of $u_i(1)$, we use the condition $\frac{\partial u}{\partial t}(x,0)=g(x)$. In deed, using $m=0$ we get $$\frac{\partial u}{\partial t}=\frac{u_i(1)-u_i(-1)}{2\Delta t} =g(x)$$ This gives us $$u_i(1)-u_i(-1)=2\Delta t g(x)$$ Also, using $m=0$ in $$u_i(m+1) =\rho^2 u_{i-1}(m)+2(1-\rho^2)u_i(m)+\rho^2 u_{i+1}(m)-u_i(m-1)$$ and rearranging few terms, we get $$u_i(1)+u_i(-1)=\rho^2 f(x_{i-1})+2(1-\rho^2)f(x_i)+\rho^2 f(x_{i-1})$$ Solve the above equations for $u_i (1)$, we get $$u_i(1)=\frac{1}{2}\rho^2(f(x_{i-1}))+(1-\rho^2)f(x_i)+\frac{1}{2}\rho^2 f(x_{i+1})+\Delta t g(x_i)$$ These equations are valid only for $\rho \leq 1$. The following Python program can be used to solve the one dimensional wave equaiton. For our application, we will use $f(x) = 2x$ if $0 < x < 1/2$ and $f(x) =2(1-x)$ if $1/2 < x < 1$. Also, we use $g(x) =0$, $n=4$, $\Delta x = 1/4$, and $\rho = 1.$

---

Two Dimensional Laplace Equaiton

Consider the two dimensional Laplace Equation with the boundary conditions

$$\frac{\partial ^2 u}{\partial x^2}+\frac{\partial ^2 u}{\partial x^2}=0$$ $$0 < x < 1, \quad 0 < y < 1$$ $$u(0,y) = 0, \quad u(x,0)=0, \quad u(x,1)=50, \quad u(1,y)=50$$ The difference equaiton becomes $$\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}=\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta y)^2}, \quad i=1, 2, \cdots, n-1$$ Let $$\Delta x = \Delta$$ $$u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}=0$$ For our purpose, let's use $n=4$. Then we have $\Delta x =\frac{1}{4}$. We will have one equation for each of the nine $(i,j)$'s as follows: $$u_{2,1}+u_{1,2}-4u_{1,1} =0$$ $$u_{3,1}+u_{1,1}+u_{2,2}-4u_{2,1}=0$$ $$u_{4,1}+u_{2,1}+u_{3,1}-4_{3,1}=0$$ $$u_{2,2}+u_{1,3}-4u_{1,2}=0$$ $$u_{3,2}+u_{1,2}+u_{2,3}+u_{2,1}-4u_{2,2}=0$$ $$u_{2,2}+u_{3,3}+u_{3,1}-4u_{3,2}=-25$$ $$u_{2,3}+u_{3,2}-4u_{1,3} = -50$$ $$u_{3,3}+u_{1,3}+u_{2,2}-4u_{2,3} = -50$$ $$u_{2,3}+u_{3,2}-4u_{3,3} = -75$$ The above system of 9 equations can be put as a matrix vector equation as $$A u =f,$$ where $$A = \begin{bmatrix} -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \\ 1 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 \end{bmatrix}$$ $$u = \begin{bmatrix} u_{1,1}\\ u_{2,1}\\ u_{3,1}\\ u_{1,2}\\ u_{2,2}\\ u_{3,2}\\ u_{1,3}\\ u_{2,3}\\ u_{3,3} \end{bmatrix}$$ $$f = \begin{bmatrix} 0\\ 0\\ 0\\ 0\\ 0\\ -25\\ -50\\ -50\\ -75 \end{bmatrix}$$ The solution is $$u = A^{-1} f$$ . The following Python program solves this matrix equation.

---

---