Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
4  
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? –  user21 Jun 19 '13 at 7: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
add comment

2 Answers 2

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.

share|improve this answer
    
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
4  
@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
    
Thank you eric for the brilliant solution –  Pawan Apr 19 at 19:42
add comment

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

share|improve this answer
    
This is nice, but is it a coupled second-order system? I think that would be more helpful to the original poster. –  Verbeia Jun 18 '13 at 23:23
    
+1 for the nice plots! –  sebhofer Jun 20 '13 at 7:45
    
Appreciate it ^^ Thank you ! –  Sektor Jun 20 '13 at 10:11
    
Thank you sector for your brilliant solution. –  Pawan Apr 19 at 19:44
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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