MATLAB:Ordinary Differential Equations

From PrattWiki
Revision as of 20:30, 10 February 2015 by DukeEgr93 (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

Introduction

This page is aimed at introducing techniques for solving initial-value problems involving ordinary differential equations using MATLAB. Specifically, it will look at systems of the form:

\( \begin{align} \frac{dy}{dt}&=f(t, y, C) \end{align} \)

where \(y\) represents an array of dependent variables, \(t\) represents the independent variable, and \(C\) represents an array of constants. Note that although the equation above is a first-order differential equation, many higher-order equations can be re-written to satisfy the form above.

In addition, the examples on this page will assume that the initial values of the variables in \(y\) are known - this is what makes these kinds of problems initial value problems (as opposed to boundary value problems).

Solving initial value problems in MATLAB may be done with two coding components. The first will be a function that accepts the independent variable, the dependent variables, and any necessary constant parameters and returns the values for the first derivatives of each of the dependent variables. In other words, you will need to write a .m function that takes \(t\), \(y\), and possibly \(C\) and returns \(f(t, y, C)\). Note that the output needs to be returned as a column vector.

The second file will be a script or function that uses the first function in concert with MATLAB's ODE solvers to calculate solutions over a specified time range assuming given initial conditions.

Formal Organization

Differential Function File

An easy way to set up the function that will calculate the values of the first derivatives of each of the dependent variables is to pass it three inputs - the independent variable, the dependent variable, and the constants. The function should return the values of the derivatives in a single matrix. Note that for multi-variable systems, the derivatives should be returned as a column vector. For many systems, this code only needs to have one (possibly algebraically complicated) command in it. A template is available at MATLAB:Ordinary Differential Equations/Templates while several examples - combined with Execution File examples - are at MATLAB:Ordinary Differential Equations/Examples.

Execution Script

Once the function for the differential is done, you need to write code to actually use it for a specific case. This file needs to set up your time base (either the initial and final time for which to solve or the complete set of times at which to solve), define the initial values for your dependent variables, and set any constants needed for the problem. You will then use one of MATLAB's ODE solvers to calculate the values for the variables at the specified times. Typically, you will then generate a graph of the answer. A template is available at MATLAB:Ordinary Differential Equations/Templates while several examples - combined with Differential File examples - are at MATLAB:Ordinary Differential Equations/Examples.

One-File Example

If the differential equation is relatively simple, the whole process can be done in a single file; for example, if you want to look at the value of some output \(y(t)\) obtained from a differential system

\( \begin{align} \frac{dy(t)}{dt}&=x(t)-\frac{y(t)}{4}\\ x(t)&=\cos(3t)\\ y(0)&=5 \end{align}\)

you can write that code as:

clear

tspan        = linspace(0, 15, 1000);
yinit        = 5;
x            = @(t) cos(3*t);
[tout, yout] = ode45(@(t, y) x(t)-y(1)/4, tspan, yinit);

plot(tout, x(tout), 'k-', tout, yout, 'k--');
legend('Input', 'Output');

Note that the \(x(t)\) is explicitly given as a function of \(t\) in the ode45 program but that the value of the variable of which you are taking a derivative is given as y(1). This is because y(1) represents the first (and only) variable of which you are solving the ODE.

As an aside, to get the same graph in Maple, you could use:

restart;

deqn := diff(y(t), t) = x(t)-(1/4)*y(t);
  
Src := x(t) = cos(3*t);

MySoln := dsolve({subs(Src, deqn), y(0) = 5}, [y(t)]);

plot(subs(Src, MySoln, [x(t), y(t)]), t = 0 .. 15);

Questions

Post your questions by editing the discussion page of this article. Edit the page, then scroll to the bottom and add a question by putting in the characters *{{Q}}, followed by your question and finally your signature (with four tildes, i.e. ~~~~). Using the {{Q}} will automatically put the page in the category of pages with questions - other editors hoping to help out can then go to that category page to see where the questions are. See the page for Template:Q for details and examples.

External Links

References