In this tutorial, we are going to discuss a MATLAB solver '

*pdepe*' that is used to solve**p**artial**d**ifferential**e**quations (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,

**is a function of sinusoidal***A**, meaning that***x***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**. Before we move forward with the coding, we need to understand first how the '***t*⩾ 0*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:Here,

*is the independent spatial variable.*

**x***is the independent time variable.*

**t***is the dependent variable being differentiated with respect to*

**y***and*

**x***. It is a two-element vector where*

**t***y(1)*is

*y₁(x,t)*and

*y(2)*is

*y₂*

*(x,t)*.

**is the symmetry constant. For cartesian coordinate, the value is**

*m**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***, and***f***refer to the coefficients in the above two PDE equations (1) and (2), which are required in a form that is usually expected by '***s***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:__

__MATLAB Program:__

**function system_of_PDEs**

**close all**

**clear**

**clc**

**% 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);**

**figure(1)**

**surf(x,t,y1)**

**title('y_1(x,t)')**

**xlabel('Distance x')**

**ylabel('Time t')**

**figure(2)**

**surf(x,t,y2)**

**title('y_2(x,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];**

**end**

**% Initial Conditions**

**function u0 = pde_ics(x)**

**u0 = [1; 0];**

**end**

**% 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];**

**end**

**end**

__Results:__

__Results:__

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.#PDE #Matlab #FiniteDifference #HyperbolicEquation #PDEPE #Blog #Blogger

Thank you for explaining the matlab pde command with examples.

ReplyDeleteThanks Rohan.

DeleteThis comment has been removed by the author.

ReplyDeleteHi Ovy,

ReplyDeleteI was working on modeling the Euler Bernoulli beam using the pdepe function in MATLAB. But in my equations, there are the second partial derivative with respect to t of the function w(Y,t), and some fourth partial derivative with respect to Y of the function w(Y,t). I searched for a similar problem, but it seems like no one has done it yet. Can you please help me with it? Thank you so much! You can take a look at the paper "Optimal Motion Planning and Energy-Based Control of a Single Mast Stacker Crane" written by Hubert Rams, Markus SchÃ¶berl, and Kurt Schlacher, equations (5a) to (5e) to visualize it easier.

Best regards.