Accelerating the pace of engineering and science

# MATLAB

## Differential Equations

This example shows how to formulate and solve three different types of differential equations using MATLAB®. MATLAB offers several numerical algorithms to solve a wide variety of differential equations.

Initial Value Problem

VANDERPOLDEMO is a function that defines the van der Pol equation.

type vanderpoldemo

function dydt = vanderpoldemo(t,y,Mu)
%VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO.

% Copyright 1984-2002 The MathWorks, Inc.

dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];


The equation is written as a system of two first order ODEs. These are evaluated for different values of the parameter Mu. For faster integration, we choose an appropriate solver based on the value of the parameter Mu.

For Mu = 1, any of the MATLAB ODE solvers can solve the van der Pol equation efficiently. The ODE45 solver used below is one such example. The equation is solved in the domain [0, 20].

tspan = [0, 20];
y0 = [2; 0];
Mu = 1;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode45(ode, tspan, y0);

% Plot of the solution
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')


For larger magnitudes of Mu, the problem becomes stiff. Special numerical methods are needed for fast integration. ODE15S, ODE23S, ODE23T, and ODE23TB can solve stiff problems efficiently.

Here is a solution to the van der Pol equation for Mu = 1000 using ODE15S.

tspan = [0, 3000];
y0 = [2; 0];
Mu = 1000;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode15s(ode, tspan, y0);

plot(t,y(:,1))
title('van der Pol Equation, \mu = 1000')
axis([0 3000 -3 3])
xlabel('t')
ylabel('solution y')


Boundary Value Problems

BVP4C solves boundary value problems for ordinary differential equations.

The example function TWOODE has a differential equation written as a system of two first order ODEs.

type twoode

function dydx = twoode(x,y)
%TWOODE  Evaluate the differential equations for TWOBVP.
%

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2002 The MathWorks, Inc.

dydx = [ y(2); -abs(y(1)) ];


TWOBC has the boundary conditions for TWOODE.

type twobc

function res = twobc(ya,yb)
%TWOBC  Evaluate the residual in the boundary conditions for TWOBVP.
%

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2002 The MathWorks, Inc.

res = [ ya(1); yb(1) + 2 ];


Prior to using BVP4C, we have to provide a guess for the solution we want represented at a mesh. The solver then adapts the mesh as it refines the solution.

BVPINIT assembles the initial guess in the form that the solver BVP4C will need. For an initial mesh of [0 1 2 3 4] and a constant guess of y(x) = 1, y'(x) = 0, call BVPINIT like this:

solinit = bvpinit([0 1 2 3 4],[1; 0]);


With this initial guess, we can solve the problem with BVP4C.

The solution sol (below) is then evaluated at points xint using DEVAL and plotted.

sol = bvp4c(@twoode, @twobc, solinit);

xint = linspace(0, 4, 50);
yint = deval(sol, xint);
plot(xint, yint(1,:),'b');
xlabel('x')
ylabel('solution y')
hold on


This particular boundary value problem has exactly two solutions. The other solution is obtained for an initial guess of

   y(x) = -1, y'(x) = 0

and plotted as before.

solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);

xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
hold off


Partial Differential Equations

PDEPE solves partial differential equations in one space variable and time.

The examples PDEX1, PDEX2, PDEX3, PDEX4, PDEX5 form a mini-tutorial on using PDEPE. Browse through these functions for more examples.

This example problem uses functions PDEX1PDE, PDEX1IC, and PDEX1BC.

PDEX1PDE defines the differential equation.

type pdex1pde

function [c,f,s] = pdex1pde(x,t,u,DuDx)
%PDEX1PDE  Evaluate the differential equations components for the PDEX1 problem.
%

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2002 The MathWorks, Inc.

c = pi^2;
f = DuDx;
s = 0;


PDEX1IC sets up the initial conditions.

type pdex1ic

function u0 = pdex1ic(x)
%PDEX1IC  Evaluate the initial conditions for the problem coded in PDEX1.
%

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2002 The MathWorks, Inc.

u0 = sin(pi*x);


PDEX1BC sets up the boundary conditions.

type pdex1bc

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
%PDEX1BC  Evaluate the boundary conditions for the problem coded in PDEX1.
%

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2002 The MathWorks, Inc.

pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;


PDEPE requires x, the spatial discretization, and t, a vector of times at which you want a snapshot of the solution. We solve this problem using a mesh of 20 nodes and request the solution at five values of t. Finally, we extract and plot the first component of the solution.

x = linspace(0,1,20);
t = [0 0.5 1 1.5 2];
sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

u1 = sol(:,:,1);

surf(x,t,u1);
xlabel('x');
ylabel('t');
zlabel('u');