Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

How to Solve a System of Partial Differential Equations in MATLAB?

In this tutorial, we are going to discuss a MATLAB solver 'pdepe' that is used to solve partial differential equations (PDEs). Let us consider the following two PDEs that may represent some physical phenomena. Sometimes, it is quite challenging to get even a numerical solution for a system of coupled nonlinear PDEs with mixed boundary conditions. This example demonstrates how we may solve a system of two PDEs simultaneously by formulating it according to the MATLAB solver format and then, plotting the results. This method is relatively easier and saves time while coding.

∂y₁/∂t = 0.375 ∂²y₁/∂x² + A(y₁ - y₂)          (1)

∂y₂/∂t = 0.299 ∂²y₂/∂x² - A(y₁ - y₂)          (2)

Where, A is a function of sinusoidal x, meaning that A = sin (x). The above two equations have two derivative terms. The time derivative, generally represents parabolic equation, and the spatial derivative defines elliptic equation. So, we have both forms of PDEs and the highest order is two in the space. There is one limitation of the 'pdepe' that you need to have the parabolic term in the PDEs in order to solve by the 'pdepe'. Now, let's assume, we have the following initial conditions:

y₁(x,0) = 1

y₂(x,0) = 0

The boundary conditions are,

∂/∂x y₁(0,t) = 0

y₂(0,t) = 0

∂/∂x y₂(1,t) = 0

y₁(1,t) = 1

The initial and boundary conditions are true if 0 ⩽ x ⩽ 1 and t ⩾ 0. Before we move forward with the coding, we need to understand first how the 'pdepe' solver accepts the PDEs in MATLAB, which form it recognizes. The general form of PDEs that the solver understands is of the following form:

The form of PDEs that the 'pdepe' solver understands


x is the independent spatial variable.

t is the independent time variable.

y is the dependent variable being differentiated with respect to x and t. It is a two-element vector where y(1) is y₁(x,t) and y(2) is y₂(x,t).

m is the symmetry constant. For cartesian coordinate, the value is 0; for cylindrical coordinate, the value is 1; and for spherical coordinate, it is 2. For our problem, it is 0 as our coordinate system is simply cartesian.

The functions c, f, and s refer to the coefficients in the above two PDE equations (1) and (2), which are required in a form that is usually expected by 'pdepe' solver.

So, in the above mentioned form, the coefficients of the PDEs may be arranged in a matrix and the equations become as,

Now, we are ready to begin the coding. In the MATLAB, we may create three functions, for example, the first function is for the equations, the second function is for the initial conditions, and the third function is for the boundary conditions. Also, we may incorporate all these functions inside another global function that is convenient as we would have everything in a single file. As the MATLAB solvers use the finite difference approach, the time integration is done with the MATLAB ‘ode15s’ solver. So, the ‘pdepe’ takes advantage of the capabilities of the stiff ‘ode15s’ solver for solving the differential-algebraic equations (DAE), which may arise when the PDEs contain elliptic equations. It is also used for handling the Jacobians with a certain sparsity pattern. Now, if there is an error during the solution process, you may try to refine the mesh size as the initial conditions are sensitive and sometimes are not consistent with the mesh size and the solver.

MATLAB Program:

function system_of_PDEs
close all

% Discretization of the simulation domain
x = linspace(0,10,100);
t = linspace(0,10,100);

% Application of the Matlab partial differential equation solver 'pdepe'
m = 0;
sol = pdepe(m,@pde_func,@pde_ics,@pde_bcs,x,t);

% The solution matrices
y1 = sol(:,:,1);
y2 = sol(:,:,2);

xlabel('Distance x')
ylabel('Time t')

xlabel('Distance x')
ylabel('Time t')

% Equations to solve
function [c,f,s] = pde_func(x,t,y,dydx) 
% Equations arranged for the 'pdepe' solver
c = [1; 1];
f = [0.375; 0.299].*dydx;
A = sin(x);
s = [A; -A];

% Initial Conditions
function u0 = pde_ics(x) 
u0 = [1; 0];

% Boundary Conditions
function [pl,ql,pr,qr] = pde_bcs(xl,ul,xr,ur,t) 
pl = [0; ul(2)];
ql = [1; 0];
pr = [ur(1)-1; 0];
qr = [0; 1];



The following results show the variations of the two dependent variables, y₁, and y₂, with respect to both space (x) and time (t). We may use the MATLAB 'Surface Plot' feature to do that.

Showing the variation of the variable y1 with respect to space and time

Showing the variation of the variable y2 with respect to space and time

#PDE #Matlab #FiniteDifference #HyperbolicEquation #PDEPE #Blog #Blogger