5
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • 5
    $\begingroup$ 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? $\endgroup$
    – user21
    Jun 19 '13 at 7:11
  • $\begingroup$ 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. $\endgroup$
    – Pawan
    Jun 19 '13 at 15:55
  • 1
    $\begingroup$ @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. $\endgroup$
    – Eric Brown
    Jun 19 '13 at 16:50
  • $\begingroup$ @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. $\endgroup$
    – Eric Brown
    Jun 19 '13 at 17:04
18
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Pawan
    Jun 19 '13 at 16:07
  • 5
    $\begingroup$ @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. $\endgroup$
    – bill s
    Jun 19 '13 at 17:21
  • 1
    $\begingroup$ Thank you eric for the brilliant solution $\endgroup$
    – Pawan
    Apr 19 '14 at 19:42
4
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • $\begingroup$ This is nice, but is it a coupled second-order system? I think that would be more helpful to the original poster. $\endgroup$
    – Verbeia
    Jun 18 '13 at 23:23
  • $\begingroup$ Appreciate it ^^ Thank you ! $\endgroup$
    – Sektor
    Jun 20 '13 at 10:11
  • $\begingroup$ Thank you sector for your brilliant solution. $\endgroup$
    – Pawan
    Apr 19 '14 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.