5

Does anybody know if Mathematica has an analogue of MATLAB's ode45 command? I need to solve a second order coupled ODE system of equations.

Aditya
  • 315
  • 1
  • 8
Pawan
  • 93
  • 1
  • 8
  • 5
    why do you need a specific integrator? The automatic method selection of Mathematica does this for you in a very sophisticated manner for a much wider range of DEs where it is not easy to know in advance which combination of integrators are needed. Only if this method selection does not work to your satisfaction, then I'd fiddle with giving specific integrators. What does NDSolve return for your specific system of ODEs? –  Jun 19 '13 at 07:11
  • NDSolve returns an interpolating function(with the domain and range) which I can only plot, but it doesn't return the result as a function. Is there any command to get the function from an interpolating function. – Pawan Jun 19 '13 at 15:55
  • 1
    @ruebenko I manually select the integrator a) for pedagogical purposes and b) to (try to) reproduce the results of someone else's integration scheme (e.g. matlab in this case). For my own stuff I let Mathematica handle it automatically. Also, sometimes I have to "port" to Fortran, where I may have LSODA at tops, sometimes less. I need to know how well the less sophisticated methods perform. – Eric Brown Jun 19 '13 at 16:50
  • @user4402 This is a numerical solution to the ODE, not an analytic solution. If you want to create your own approximate representation, then you can extract the data from the InterpolatingFunction itself, or use it to generate the data at points that you need. – Eric Brown Jun 19 '13 at 17:04

2 Answers2

18

According to the Mathematica documentation on this page:

Here is how to define a 5(4) pair of Dormand and Prince coefficients [DP80]. This is currently the method used by ode45 in MATLAB.

DOPRIamat = {
   {1/5},
   {3/40, 9/40},
   {44/45, -56/15, 32/9},
   {19372/6561, -25360/2187, 64448/6561, -212/729},
   {9017/3168, -355/33, 46732/5247, 49/176, -5103/18656},
   {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84}};
DOPRIbvec = {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0};
DOPRIcvec = {1/5, 3/10, 4/5, 8/9, 1, 1};
DOPRIevec = {71/57600, 0, -71/16695, 71/1920, -17253/339200, 
   22/525, -1/40};
DOPRICoefficients[5, p_] :=
  N[{DOPRIamat, DOPRIbvec, DOPRIcvec, DOPRIevec}, p];

Then:

NDSolve[system, 
        Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5, 
                   "Coefficients" -> DOPRICoefficients, "StiffnessTest" -> False}]

where

system

is the second order ODE specified in the usual Mathematica manner for ODE's.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Eric Brown
  • 4,406
  • 1
  • 18
  • 36
  • The simple NDSolve solves the ODE's i am using but it returns an interpolationg function. Whereas I want the function itself so that I can make further changes to it. How should I get the function from the interpolating function. – Pawan Jun 19 '13 at 16:07
  • 5
    @user4402 -- Matlab's ODE45 returns a vector of values and the times at which those values occur. If you wish to duplicate this output, you can simply evaluate the interpolating function at the points you wish. For example, with fInt the interpolating function, Table[fInt[x],{x,0,10,0.1}] will return the values at the points 0, 0.1, 0.2, ... 10. – bill s Jun 19 '13 at 17:21
  • 1
    Thank you eric for the brilliant solution – Pawan Apr 19 '14 at 19:42
4

Check the documentation centre right here: LINK

And a simple example (modelling coupled springs with added nonlinear restoring forces (for large vibrations)):

eqn = {x''[t] == 
0.4 x[t] + -1/6 x[t]^3 - 1.808 (x[t] - y[t]) - 
 1/10 (x[t] - y[t])^3, 
y''[t] == -1.808 (y[t] - x[t]) - 1/10 (y[t] - x[t])^3, 
x[0] == -0.6, x'[0] == 1/2, y[0] == 3.001, y'[0] == 5.9};

sol = NDSolve[eqn, {x, y}, {t, 0, 200}]

ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 200}]

x[t] vs y[t] Plot

ParametricPlot[Evaluate[{x[t], x'[t]} /. sol], {t, 0, 200}]

x[t] vs x'[t] Phase Plot

ParametricPlot[Evaluate[{y[t], y'[t]} /. sol], {t, 0, 200}]

y[t] vs y'[t] Phase Plot

Sektor
  • 3,320
  • 7
  • 27
  • 36