The Numerical Solution of Partial Differential Equations with an Example. C Code and Octave Script
1. Abbreviation
2. Introduction
Most partial differential equations cannot be solved analytically (in exact) form, therefore numerical (approximate) solution methods are used. In this article I want to consider an example of a numerical solution of a simple partial differential equation of two independent variables.
In general a second-order linear partial differential equation can be written:
Equation (10) is called parabolic in the case of:
Equation (10) is called hyperbolic in the case of:
Equation (10) is called elliptical in the case of:
As an example we will consider the hyperbolic wave equation describing transverse vibrations of a string, longitudinal vibrations of a rod, electrical vibrations in a wire and etc.:
Let the desired function u(x,t) be given in the area:
To obtain an unique solution, initial and boundary conditions are set:
This task will be solved for zero boundary conditions:
The numerical solution method is based that continuous independent variables are replaced by a grid of discrete values, and partial derivatives are replaced by difference equations defined on this grid. Not every choice of a grid of discrete values and difference equations leads to a stable solution of the partial differential equation. An accumulating calculation error can lead to an incorrect result. Therefore it is necessary to first investigate grid model, difference equations and prove the stability of its solution. See the [1], [2]
Note:
This wave equation (14) – (17) can also be solved analytically. Having an accurate (analytical) solution, it will be possible to estimate the accuracy of the numerical solution.
3. The Analytical Solution of the Wave Partial Differential Equation
To solve equation (14) – (17), the method of separation of variables is used. See also the [3]. A partial solution of the equation will be looked for in the form:
Then:
Substitute (19) into (14):
Note: it can be proved that for θ ≤ 0 we get a trivial (null) solution of no interest.
The “θ” does not depend on the variables “x”, “t”. As a result we obtain two ordinary linear differential equations:
Solving (22):
where “A”, “B”, “C”, “D” are constants that will be defined from (16) and (17a)
Using (17a) and (23):
Then from (24):
A = 0 (26)
From (25):
Then:
Rewrite (23) taking into account (26), (27):
Note
In this case the constant “B” is included in “C” and “D”
The general solution of equation (14) can be written as the sum of the simple solutions (28):
The unknown “C” and “D” are calculated from the initial conditions (16):
Thus these coefficients are found by calculating the coefficients of the Fourier series:
Having a general solution to the wave equation, consider a specific example:
Initial conditions:
Boundary conditions:
Analytical solution for (34), (35a), (35b), (36):
Figure 3-1 shows the solution (37):
Figure 3-1: Analytical Solution of the (34) – (37)
4. The Numerical Solution of the Wave Partial Differential Equation
The desired function u(x,t) is replaced by a function defined on a grid of discrete values:
where
h – step along the x-axis
Ts – step along the t-axis
Consider the following explicit grid stencil of the finite difference method for the wave equation (14) (see also [1]):
Figure 4-1: The explicit method stencil
Then the wave equation (14) can be written using finite differences / difference equation:
This method is called explicit because the unknown value from the next layer “n” is calculated using the values already obtained from the previous layer. The accuracy of the approximation of the second partial derivatives is O(h^2), O(Ts^2). See the Appendix.
Let’s determine the step values of the “h” and “Ts” at which calculations using the formula (40) give stable solutions. To do this, we will use the Fourier stability analysis using the von Neumann condition. The method is valid for linear partial differential equations that are approximated by finite differences. See also [2].
For analysis a harmonic function is taken, which is the solution of a linear difference equation (an eigenfunction of a linear difference operator) of the form:
The difference equation will be stable if for any frequency ω the value |λ| ≤ 1 (43). This criterion is called the von Neumann condition. This means that any disturbances caused by errors in rounding, approximation or initial conditions will decrease when calculating subsequent layers.
Substitute (41) in (39):
Let’s introduce the notation:
Then:
Let’s find the roots of this equation (45):
From Vieta’s formulas:
Consider the value:
In this case the solution of the equation is two complex conjugate numbers:
Thus the von Neumann conditions (43) are fulfilled and the difference equation (40) gives a stable solution.
When the value is:
The roots of the equation can take values greater than one in absolute magnitude, which leads to an unstable solution of the difference equation. For example:
Then:
In this case one root takes the value:
The von Neumann condition is not fulfilled, which leads to an unstable solution (40) under condition (51).
Thus we proved that the difference equation (40) gives a stable solution under the condition (49). Now we solve the problem (34) – (36) numerically using (40). To do this, select the steps:
h = 0.5m
Ts = 0.1s (56)
Let’s check the condition (49):
Condition (49) is fulfilled.
To start the formula (40), we need to know the first two layers of the u(kh, 0) and u(kh, Ts). The zero layer is known from the initial conditions (35a):
The next first layer is calculated using Truncation of the Taylor series (2) and initial conditions (35b):
The accuracy of the calculation is O(Ts^2) which follows from the remainder term of the Taylor series. At the same time the accuracy of calculations of the first layer can be improved if the second derivative is known:
Then (59) can be rewritten:
The accuracy of calculations will improve in this case to O(Ts^3).
Next using the boundary conditions (36) and the formula (40), the remaining layers are calculated:
Figure 4-2 shows the numerical solution (34)-(36):
Figure 4-2: Numerical Solution of the (34) – (36)
See the NumericalPartialDifferentialEquation.m and NumericalPartialDifferentialEquation.c/h
5. Conclusions
-
For each partial differential equation, it is necessary to build a finite difference stencil and the difference equation of which has a stable solution with a certain choice of grid steps and provides the necessary accuracy of the solution
-
The numerical example above gives a maximum calculation error of 2%. See the NumericalPartialDifferentialEquation.m
If higher accuracy is required, then it is necessary to reduce the steps “h”, “Ts”
-
Fourier stability analysis is valid for linear difference equations. If the equation is nonlinear, then it is necessary to linearize the equation and then use Fourier stability analysis
-
In our above example, an explicit stencil of the difference equation was used. Here is an example of an implicit stencil:
Figure 5-1: The implicit method stencil
Then the difference equation will be written as:
In (63) the 3 values from layer “n” are unknown at once:
To find these values (65), linear algebraic equations are compiled, which are then solved by numerical methods. See the [1], [2]
It can be proved that when:
The difference equation (63) gives a stable solution for any choice of steps “h” and “Ts”, while the accuracy corresponds to O(h^2), O(Ts^2).
6. Octave GNU file PartialNumericalDifferentialEquation.m
The .m script generates the analytical and numerical solutions of the partial differential equation (34)-(36) and result plots. To do this, the following functions are defined in this file:
% Numerical Solution Of the Partial Differential Equation with the Explicit Method
% u_tt(x,t) = a_speed^2*u_xx(x,t)
% Input parameters:
% u_0_t – boundary values
% u_L_t – boundary values
% u_x_0 – initial values
% u_t_x_0 – initial values
% u_tt_x_0 – initial values
% a_speed – Partial Differential Equation Parameter
% h – sample step
% Ts – sample time
% L – x space: 0…L with step Ts
% T – time space: 0…T with step h
% Output parameter:
% u_output – output
function [u_output,u_pattern] = NumericalSolutionOfPartialDifferentialEquation(u_0_t, u_L_t, u_x_0, u_t_x_0, u_tt_x_0, a_speed, h, Ts, L, T)
% Float Table to File
% Support C-Code
% Input parameters:
% ParameterVector – data array
% FileNameString – output file name
function FloatParamVector2file(ParameterVector, FileNameString)
7. C Language PartialNumericalDifferentialEquation.c/h
The demo SW is realized the numerical solutions of the partial differential equation (34)-(36). See the files PartialNumericalDifferentialEquation.c/h
Functions:
/*
Wave Equation: u_tt(x,t) = a_speed^2*u_xx(x,t)
Partial Differential Equation Initialization
Initialization of the working arrays
Input: array_ptr – 2-Dimensional array pointer
Return: void
*/
void PartialDifferentialEquationInit(float* array_ptr)
/*
Wave Equation: u_tt(x,t) = a_speed^2*u_xx(x,t)
Numerical Solution Of Partial Differential Equation
Input
Input: array_ptr – 2-Dimensional array pinter
Return: void
*/
void NumericalSolutionOfPartialDifferentialEquation(float* array_ptr)
8. Download the PartialNumericalDifferentialEquation.c/h/m
You can download the files:
NumericalPartialDifferentialEquation.m
NumericalPartialDifferentialEquation.c
NumericalPartialDifferentialEquation.h
with the button:
9. Literature / References
[1] Kireev V.I., Panteleev A.V. “Numerical methods in examples and tasks”, Moscow, Publishing House Higher School, 2008, ISBN 978-5-06-004763-9
[2] Pirumov U.G. “Numerical methods”, Moscow, Publishing House of the Moscow Aviation Institute, 1998, ISBN 5-7035-2190-4
[3] N.S. Piskunov “Differential and integral calculus for higher technical educational institutions”, volumes 1 and 2, Moscow, “Nauka” Publishing House, 1978
Appendix 1. The accuracy of calculations of the second partial derivatives in the difference approximation (39)
Let’s estimate the accuracy of the formula (39):
To do this, the Taylor series are written with a remainder term for (k+1) h and (k-1)h:
Note:
Summation of (A2) and (A3):
Then:
