Matlab Numeric ODE Solutions
A very useful website: Mathworks online documentation
Numeric ODE Solutions
Matlab has a suite of functions to help solve ordinary differential equations using numeric techniques. The examples below are only a brief introduction to the subject.
Look up ode23 in MATLAB help. We will not use this command for our example, but this key word search will take you to the documentation for all of the ODE solver functions.
General steps for solving ODEs
- MATLAB can only solve ODEs of the first order. However, MATLAB can solve a system of first order equations. Thus, if you have an ODE of the n-th degree you must break your equation into n first order ODEs.
Define your problem in either a script file or a function file:
- Range of values to solve over. A time range, for example.
- Define the initial conditions for your first order differential equations.
- Any options you want to change (see MATLAB help). The odeset function can be used to specify options for the various ode functions
- Set up a function that takes as input a value of an independent variable and a dependent variable (t and y, for example), and returns dy/dt. If you are solving a system of equations, the y and the return value dy/dt become lists instead of single values.
- The outputs of the ode45 function are a list (column) of t values and y is a matrix of numbers that are values on the curve(s) that are the approximate solution(s) to the ODE equation(s). y will have as many columns as there are first order equations being solved.
- Solve dy/dt = cos(t) with ode45 for y(0)=2 and the time range 0 to 2π. Plot the analytic and numeric solution on the same plot. Examine the differences. The analytic solution is y=sin(t)+2.
- Solve a 2nd order equation x'' + 16x = sin(4.3t), x(0)=0 and x'(0)=0, time range 0 to 2π.
To do this you must re-write the equation as a system of two first order equations. Use y=x and z=x' to get z'=sin(4.3t)-16y and y'=z. You must create a separate function with a given format for use by ode45.
- The equations of motion you learned in EF 151 are simple cases of a solution to an ODE. Use y''=-g and ode45 to plot the height of a projectile with respect to time for y(0)=42 and y'(0)=12. If you add air resistance the equation becomes y''=-g ± c*y'2 (drag is proportional to velocity squared, drag direction is opposite velocity). Add a second curve to the plot that incorporates air resistance with c=0.01. Modify to terminate the solution when y=0.
- Extra - modify the solution above for 2-D motion.
- Railgun Example
Foucault Pendulum Example
Here is another example of the use of ODE45.File(s):
ReferencesMatlab Demystified, David McMahon, Chapters 7 and 9
Numerical Computing with Matlab, Cleve Moler