Where μ = 1200, t is from 0 to 7000, and ODE15s is to be used.
tfinal = 7000;
delta_t = 0.01;
t = [0:delta_t:tfinal];
initial = [1 1]; % initial conditions
[t,w] = ode15s(@sub2,t,initial); % line with error
figure
plot(t,w,'LineWidth',2)
%%%%%%%%%% Subfunction
function d_dt = sub2(t,w)
mu = 1200;
% turning d^2y/dt^2 = mu(1-y^2)dydt - y into dzdt = mu(1-y)^2 - y
% dydt = z
z = w(1);
y = w(2);
d_dt = zeros(2, 1); % initialize the column vector of equations
d_dt(1,1) = mu*z*(1-y)^2 - y; % first equation
d_dt(2,1) = z; % second equation
end
I made the code above, but i keep getting the error:
Warning: Failure at t=1.702777e-01. Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (4.440892e-16) at time t.
> In ode15s (line 730)
In sec_ord_diff (line 5)
I think it has something to do with μ being too large, or maybe the time step is too small/ big. Can anyone help me fix the issue? Thank you

 채택된 답변

Star Strider
Star Strider 2019년 10월 27일

0 개 추천

Your ‘sub2’ function is not correct.
sub2 = @(t,w,mu) [w(2); w(1)-mu*(w(1)^2-1)*w(2)];
mu = 1200;
tfinal = 7000;
delta_t = 0.1;
t = [0:delta_t:tfinal];
initial = [1 1]; % initial conditions
[t,w] = ode15s(@(t,w)sub2(t,w,mu),t,initial); % line with error
figure
semilogy(t,w,'LineWidth',2)

댓글 수: 2

Tarek Chahattou
Tarek Chahattou 2019년 10월 28일
Thank you so much! Also, the semilogy was a nice thing i didn't know about.
How did you get the w(1)-mu*(w(1)^2-1)*w(2)? When i try to simplify of of the equations, I keep getting mu*(1-w(1)^2)*w(2)-w(1). It seems to be a sign difference or something. I keep wanting to do:
dy/dt = z
dz/dt - μ(1-y^2)z + y = 0
then
dy/dt = z
dz/dt = μ(1-y^2)z - y
As always, my pleasure!
Actually, I let the Symbolic Math Toolbox do the work:
syms t y(t) mu
Eq = diff(y,2) == mu*(1-y^2)*diff(y)+y;
[VF,Sbs] = odeToVectorField(Eq);
I then couid have used the matlabFunction function, however it was easier to just transcribe the ‘VF’ results.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

질문:

2019년 10월 27일

댓글:

2019년 10월 28일

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by