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.
2 Answers
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.
-
$\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$– PawanJun 19, 2013 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 sJun 19, 2013 at 17:21 -
1
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}]
ParametricPlot[Evaluate[{x[t], x'[t]} /. sol], {t, 0, 200}]
ParametricPlot[Evaluate[{y[t], y'[t]} /. sol], {t, 0, 200}]
NDSolve
return for your specific system of ODEs? $\endgroup$