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.