Differential Equations

Open Live Script

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

Initial value problems
Boundary value problems
Delay differential equations
Partial differential equations

Initial Value Problem

vanderpoldemois a function that defines the van der Pol equation

$\frac{d^{2} y}{d t^{2}} - μ (1 - y^{2}) \frac{d y}{d t} + y = 0 .$

typevanderpoldemo

function dydt = vanderpoldemo(t,y,Mu) %VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO. % Copyright 1984-2014 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 ordinary differential equations (ODEs). These equations are evaluated for different values of the parameter $μ$ . For faster integration, you should choose an appropriate solver based on the value of $μ$ .

For $μ = 1$ , any of the MATLAB ODE solvers can solve the van der Pol equation efficiently. Theode45solver is one such example. The equation is solved in the domain $[0, 20]$ with the initial conditions $y (0) = 2$ and $\frac{dy}{dt} |_{t = 0} = 0$ .

tspan = [0 20]; y0 = [2; 0]; Mu = 1; ode = @(t,y) vanderpoldemo(t,y,Mu); [t,y] = ode45(ode, tspan, y0);% Plot solutionplot(t,y(:,1)) xlabel('t') ylabel('solution y') title('van der Pol Equation, \mu = 1')

Figure contains an axes object. The axes object with title v a n blank d e r blank P o l blank E q u a t i o n , blank mu blank = blank 1 contains an object of type line.

For larger magnitudes of $μ$ , the problem becomesstiff. This label is for problems that resist attempts to be evaluated with ordinary techniques. Instead, special numerical methods are needed for fast integration. Theode15s,ode23s,ode23t, andode23tbfunctions can solve stiff problems efficiently.

This solution to the van der Pol equation for $μ = 1000$ usesode15swith the same initial conditions. You need to stretch out the time span drastically to $[0, 3000]$ to be able to see the periodic movement of the solution.

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')

Figure contains an axes object. The axes object with title v a n blank d e r blank P o l blank E q u a t i o n , blank mu blank = blank 1 0 0 0 contains an object of type line.

Boundary Value Problems

bvp4candbvp5csolve boundary value problems for ordinary differential equations.

The example functiontwoodehas a differential equation written as a system of two first-order ODEs. The differential equation is

$\frac{d^{2} y}{d t^{2}} + | y | = 0 .$

typetwoode

function dydx = twoode(x,y) %TWOODE Evaluate the differential equations for TWOBVP. % % See also TWOBC, TWOBVP. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. dydx = [ y(2); -abs(y(1)) ];

The functiontwobchas the boundary conditions for the problem: $y (0) = 0$ and $y (4) = - 2$ .

typetwobc

function res = twobc(ya,yb) %TWOBC Evaluate the residual in the boundary conditions for TWOBVP. % % See also TWOODE, TWOBVP. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. res = [ ya(1); yb(1) + 2 ];

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

Thebvpinitfunction assembles the initial guess in a form you can pass to the solverbvp4c. For a mesh of[0 1 2 3 4]and constant guesses of $y (x) = 1$ and $y' (x) = 0$ , the call tobvpinitis:

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

With this initial guess, you can solve the problem withbvp4c. Evaluate the solution returned bybvp4cat some points usingdeval, and then plot the resulting values.

sol = bvp4c(@twoode, @twobc, solinit); xint = linspace(0, 4, 50); yint = deval(sol, xint); plot(xint, yint(1,:)); xlabel('x') ylabel('solution y') holdon

Figure contains an axes object. The axes object contains an object of type line.

This particular boundary value problem has exactly two solutions. You can obtain the other solution by changing the initial guesses to $y (x) = - 1$ and $y' (x) = 0$ .

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,:)); legend('Solution 1','Solution 2') holdoff

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Solution 1, Solution 2.

Delay Differential Equations

dde23,ddesd, andddensdsolve delay differential equations with various delays. The examplesddex1,ddex2,ddex3,ddex4, andddex5form a mini tutorial on using these solvers.

Theddex1example shows how to solve the system of differential equations

$\begin{array}{l} y_{1}^{'} (t) = y_{1} (t - 1) \\ y_{2}^{'} (t) = y_{1} (t - 1) + y_{2} (t - 0.2) \\ y_{3}^{'} (t) = y_{2} (t) . \end{array}$

You can represent these equations with the anonymous function

ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];

The history of the problem (for $t \leq 0$ ) is constant:

$\begin{array}{l} y_{1} (t) = 1 \\ y_{2} (t) = 1 \\ y_{3} (t) = 1 . \end{array}$

You can represent the history as a vector of ones.

ddex1hist = 1 (3,1);

A two-element vector represents the delays in the system of equations.

lags = [1 0.2];

Pass the function, delays, solution history, and interval of integration $[0, 5]$ to the solver as inputs. The solver produces a continuous solution over the whole interval of integration that is suitable for plotting.

sol = dde23(ddex1fun, lags, ddex1hist, [0 5]); plot(sol.x,sol.y); title({'An example of Wille and Baker','DDE with Constant Delays'}); xlabel('time t'); ylabel('solution y'); legend('y_1','y_2','y_3','Location','NorthWest');

Figure contains an axes object. The axes object with title An example of Wille and Baker DDE with Constant Delays contains 3 objects of type line. These objects represent y_1, y_2, y_3.

Partial Differential Equations

pdepesolves partial differential equations in one space variable and time. The examplespdex1,pdex2,pdex3,pdex4, andpdex5form a mini tutorial on usingpdepe.

This example problem uses the functionspdex1pde,pdex1ic, andpdex1bc.

pdex1pdedefines the differential equation

$π^{2} \frac{\partial u}{\partial t} = \frac{\partial}{\partial x} (\frac{\partial u}{\partial x}) .$

typepdex1pde

函数[c、f、s] = pdex1pde (x, t, u, DuDx) % pdex1pdeEvaluate the differential equations components for the PDEX1 problem. % % See also PDEPE, PDEX1. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. c = pi^2; f = DuDx; s = 0;

pdex1icsets up the initial condition

$u (x, 0) = \sin π x .$

typepdex1ic

function u0 = pdex1ic(x) %PDEX1IC Evaluate the initial conditions for the problem coded in PDEX1. % % See also PDEPE, PDEX1. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. u0 = sin(pi*x);

pdex1bcsets up the boundary conditions

$u (0, t) = 0,$

$π e^{- t} + \frac{\partial}{\partial x} u (1, t) = 0 .$

typepdex1bc

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) %PDEX1BC Evaluate the boundary conditions for the problem coded in PDEX1. % % See also PDEPE, PDEX1. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. pl = ul; ql = 0; pr = pi * exp(-t); qr = 1;

pdeperequires the spatial discretizationxand a vector of timest(at which you want a snapshot of the solution). Solve the problem using a mesh of 20 nodes and request the solution at five values oft. 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');

Figure contains an axes object. The axes object contains an object of type surface.

Differential Equations

Initial Value Problem

Boundary Value Problems

Delay Differential Equations

Partial Differential Equations

See Also

Related Topics