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.
4 个评论
Ameer Hamza
2020-6-13
Can you show the equations of your ODE in mathematical form?
Louis De Jager
2020-6-13
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
2020-6-13
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
2020-6-13
编辑：Louis De Jager
2020-6-13
采纳的回答
Star Strider
2020-6-13
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.
.
4 个评论
Louis De Jager
2020-6-13
Star Strider
2020-6-13
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
2020-6-13
Star Strider
2020-6-13
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
更多回答（1 个）
强 陈
2024-4-7
0 个投票
Hello，I am also learning Arnold's tongue recently, can I study your ODEvcnt code？
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
