Divergent solution using ode45

8 views (last 30 days)
Abhik Saha
Abhik Saha on 9 Sep 2022
Commented: Jano on 14 Oct 2022
I have a code (see below) that solve the second order differential equation. As a output it gives the position and velocity as a function of time. For some values of omega 0.48,0.5,0.85 it gives the diverging solution of position. I guess this is due to some error but I can not find the error. How to overcome this errors for above omegas. Please help regarding this. I am new to MATLAB.
%% MAIN PROGRAM %%
t=linspace(0,1000,5000);
y0=[1 0];
omega=0.85;
[tsol, ysol]=ode45(@(t,y0) firstodefun4(t,y0,omega), t, y0, omega);
position=ysol(:,1);
velocity=ysol(:,2);
%% FUNCTION DEFINITION%%
function dy=firstodefun4(t,y0,omega)
F=1;gamma=0.01;omega0=1;
kappa=0.15;
dy=zeros(2,1);
dy(1)=y0(2);
dy(2)=2*F*sin(omega*t)-2*gamma*y0(2)-omega0^2*(1+4*kappa*sin(2*omega*t))*y0(1);
end

Answers (2)

Torsten
Torsten on 9 Sep 2022
Edited: Torsten on 9 Sep 2022
We don't know the reason - technically, your code is correct. Although I would change it as shown below.
What makes you think that the solution you get is incorrect ? Do you have any indication that it should be bounded as t -> 00 ? For me, the behaviour looks very regular.
%% MAIN PROGRAM %%
t=linspace(0,50,5000);
y0=[1 0];
omega=0.85;
F=1;
gamma=0.01;
omega0=1;
kappa=0.15;
[tsol, ysol]=ode45(@(t,y0) firstodefun4(t,y0,omega,F,gamma,omega0,kappa), t, y0);
position=ysol(:,1);
velocity=ysol(:,2);
plot(tsol,position)
%% FUNCTION DEFINITION%%
function dy=firstodefun4(t,y0,omega,F,gamma,omega0,kappa)
dy=zeros(2,1);
dy(1)=y0(2);
dy(2)=2*F*sin(omega*t)-2*gamma*y0(2)-omega0^2*(1+4*kappa*sin(2*omega*t))*y0(1);
end
  4 Comments
Jano
Jano on 14 Oct 2022
Is it possible for this to accept a range of value for omega, for example omega = [0.85:5].

Sign in to comment.


Sam Chak
Sam Chak on 12 Sep 2022
If you reduce the value for , then the system should be stable.
For more info, please check out Floquet theory.
% MAIN PROGRAM %
n = 400;
t = linspace(0, n*100, n*10000+1);
y0 = [1 0];
omega = 0.85;
F = 1;
gamma = 0.01;
omega0 = 1;
kappa = 0.1469;
[tsol, ysol] = ode45(@(t, y) firstodefun4(t, y, omega, F, gamma, omega0, kappa), t, y0);
position = ysol(:, 1);
velocity = ysol(:, 2);
plot(tsol, position)
% ODE
function dy = firstodefun4(t, y, omega, F, gamma, omega0, kappa)
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = 2*F*sin(omega*t) - 2*gamma*y(2) - omega0^2*(1 + 4*kappa*sin(2*omega*t))*y(1);
end
  3 Comments
Abhik Saha
Abhik Saha on 14 Sep 2022
@Sam Chak Thank you for sharing the information related Floquet theorem. I will watch the videos and if I find some criteria then I will post here.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!