ODE45 solver, with changing initial conditions

I'm trying to numerically find the transition curves for a ODE, my code is supposed to do this by finding the solution to the ode, determining at which point the solution "blows up" and then storing the values for v and epsilon (epp) within an array.

However when running my code I keep on getting the following errors:

Unrecognized function or variable 'ODEvcnt'.

Error in ode2>@(t,y)dtheta(t,y,ODEvcnt,ODEeppcnt) (line 33)

sol = ode45(@(t,y) dtheta(t,y,ODEvcnt,ODEeppcnt),tspan,y0);

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in ode2 (line 33)

sol = ode45(@(t,y) dtheta(t,y,ODEvcnt,ODEeppcnt),tspan,y0);

I have attatched my code bellow:

Any help or advise would be much appreciated, thank you.

[Ameer Hamza]

Can you show the equations of your ODE in mathematical form?

[Louis De Jager]

The equation is as follows:

I have also attatched the analytical result of that which I'm trying to numerically determine:

Note phi and psi are the two different cases for v (being 1 and 4 respectively)

[Ameer Hamza]

Where is phi and psi in this equation? What does this graph represent? The ODE is between tau and theta, so how do you get this graph between phi and epsilon.

[Louis De Jager]

Phi and psi are two different v's with phi being for the region of v=1 and psi being fot the region of v =4 (the two cases v1 and v2)

The ode has various states: a stable state in which it produces a bounded oscillation and an unstable state during which the oscillations (of theta) grows exponentially.

Hence the graph shows the Arnold tongues for the system.

The program tries to identify the specific epsilon value for a given velocity at which the solution becomes unstable (blows up).

Answer 1

[Star Strider]

This runs without error. You need to determine if it produces the desired result:

syms y x epp t v

%cases
vM = 1; %Max's case
vMM =4; %Miu-Miu's case


%initial condition
initialY = pi/4;

%define expansions
v1 = 0.5:0.001:1.5; %max v
v2 = 3.5:0.001:4.5; % Miu Miu v
tspan = linspace(0,0.01,100); %time scale
epp = 0:0.0001:1; %epsilon possibilities

%storage variables
CritEpp = zeros(1000,2); % critical epsilon array
A = zeros(1,1000); % array to store applitudes

% sol=[t,y];                                        % 'sol' Is Not Defined At This Point

%definitiion od dtheta


for vcnt = 1:1:1000
    
    for eppcnt = 1:1:10000
        
        
        y0=[0,initialY,vcnt,eppcnt];
        ODEvcnt = vcnt;
        ODEeppcnt = eppcnt;
        
        sol = ode45(@(t,y) dtheta(t,y,ODEvcnt,ODEeppcnt),tspan,y0);
        
        A = max([sol.x.' sol.y.'],[],2);
    end
    
    if A(vcnt) > epp(eppcnt)
        CritEpp(vcnt,:) = [v1(vcnt) epp(eppcnt)];
    else
        continue;
    end
    
end


function dy = dtheta(t,y,ODEvcnt,ODEeppcnt)
    vArr = 0.5:0.001:1.5;
    eppArr = 0:0.0001:1;
    
    
    dy = -(vArr(ODEvcnt))*y - (eppArr(ODEeppcnt))*cos(2*t);
end

These were not defined previously, anywhere in your code:

        ODEvcnt = vcnt;
        ODEeppcnt = eppcnt;

I took a guess as to what they should be, based on how you use them in your code. Make the appropriate corrections if I guessed wrong.

.

[Louis De Jager]

Thank you very much.

Your guess was perfect, sorry about the confussion; would you mind just explaining what was wrong with the ODE45 function?

Was it simply because I forgot to define ODEvcnt and ODEeppcnt?

[Star Strider]

My pleasure!

‘Was it simply because I forgot to define ODEvcnt and ODEeppcnt?’

Yes. That is the reason the ode45 call failed. You are referencing the elements of the vectors in ‘dtheta’, so you need to pass the indices to it.

There were also a few other problems, specifically with respect to the ‘A’ calculation, and the difference between returning a ‘sol’ structure and the time vector and the integrated differential equations. Also, ‘CritEpp’ had the wrong dimensions.

Your essential code design is correct. You simply forgot to define and pass the vector subscript references.

[Louis De Jager]

I see, thank you very much for the help.

Enjoy the rest of your day

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

Answer 2

[强 陈]

Hello，I am also learning Arnold's tongue recently, can I study your ODEvcnt code？