# Searching for periodic solutions of Mathieu equation using MATLAB's ode45 and a crude shooting algorithm

I am numerically simulating the Mathieu equation using ODE45 and I have to keep changing the parameters delta and epsilon for each simulation to get the required peiodic solution.

Following is the MATLAB code I am running:

function xdot=trial2(t,x)
delta=0.1045;epsilon=0.0048685;
xdot=[x(2);(-delta-epsilon*cos(t))*x(1)-0.7*delta*abs(x(1))];

[t,x]=ode45('trial2',[0 10000000],[0;1]);
hold on; plot(t,x(:,1),'r');
clear all;


It is taking me 10 minutes for each simulation if I simulate for 10000000 time and it doubles if I increase the time 5 fold. After each simulation, I analyze the result and change the parameters accordingly so that I can come closer to a periodic solution. I have to do this for maybe about 30-40 parameters maybe more. I am running this on a 64 bit desktop with 8gb ram. The program runs out of memory on a 32 bit desktop.

1) Is there somehow I can improve my computational time? I might have to go from 10000000 to 50000000, maybe even more. 2) Can you help me out in implementing the parfor loop so that I can run simultaneous simulations for more number of parameters. That might be of some help too.3) WIll it help if I run this on a faster processor than i5-2400 CPU?

Thank you.

• Could you say more about your ultimate objective? Have you seen this example mathworks.com/help/matlab/ref/bvp4c.html where a Mathieu equation is solved as a boundary value problem? As you've observed, solving this equation by repeatedly solving the initial value problem for different values of the parameter is not very efficient. – Bill Greene Jul 19 '14 at 13:08

How accurate must your propagation be? You can play with the accuracy using the odeset function (to create a set of options for ode45). Since ode45 is variable step, the accuracy has a huge influence on the computational speed. The higher the required accuracy, the more steps ode45 will take. By default, the relative accuracy is set to 1e-3 and the absolute accuracy to 1e-6.

You could also consider integrating using another language (such as C or Fortran) and Mex functions or loading libraries (see loadlibrary). This will definitely result in huge speedups, since compiled languages are much faster than Matlab for this type of operations.

1) Unless the function you provided in your question changes, the period of the function should be constant. In that case, I suggest using the ability of ode45 to locate "events" -- such as zero crossings (see the documentation for odeset). That way you can simulate for a shorter stretch of time, and use the accurate estimations of zero-crossings to measure the period. That should reduce your computational time by reducing the number of solver steps required.

2) You can certainly use a parfor loop to explore parameters:

vfDelta = linspace(0, 1, 100);
vfEpsilon = linspace(0, 1, 100);
[mfD, mfE] = ndgrid(vfDelta, vfEpsilon);
mfParams = [mfD(:) mfE(:)];

parfor nParam = 1:size(mfParams, 1)
% - Set up an anonymous function to pass to ode45, that contains your parameters
fhSystem = @(t,x)trial2_params(t, x, mfParams(nParam, 1), mfParams(nParam, 2));

% - Solve the ODE
[t,x]=ode45(trial2, [0 1000], [0;1]);

% - Do your analysis here...
vfPeriod(nParam) = ...
end

function trial2_params(t, x, fDelta, fEpsilon)
...
end


3) A faster CPU will always help...

• Alternatively, you can write a small function that uses ode45 to solve the system for a given set of parameters and return the period, and then use an optimisation function like fminunc or so on to find your optimal parameter set. – Dylan Richard Muir Aug 6 '14 at 19:38