Divergent solution using ode45

13 Ansichten (letzte 30 Tage)
Abhik Saha
Abhik Saha am 9 Sep. 2022
Kommentiert: Jano am 14 Okt. 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

Antworten (2)

Torsten
Torsten am 9 Sep. 2022
Bearbeitet: Torsten am 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 Kommentare
Abhik Saha
Abhik Saha am 9 Sep. 2022
I have checked with RelTol and Abstol but it does not change much I have also check with smaller time steps but again the result remains same. If possible can you try with RADAU5 from Hairer-Wanner because I do not know how to use this. HJust to check whether it gives divergent result or not. Apart from that I have not tried with kappa=0. But I have tried with other ode integrator it does not change much.
Jano
Jano am 14 Okt. 2022
Is it possible for this to accept a range of value for omega, for example omega = [0.85:5].

Melden Sie sich an, um zu kommentieren.


Sam Chak
Sam Chak am 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 Kommentare
Sam Chak
Sam Chak am 13 Sep. 2022
@Abhik Saha, I changed the value of κ (via trial-and-error) because I thought of reducing the effect of the periodic function in the term. The original should be restored. You can definitely run the simulations for different values of angular frequency ω. There is nothing wrong about the response being divergent as (unstable).
You need to work out the algebra in order to find out fundamental matrix solution of the system, so that you can determine the stability codition following from the Floquet theorem. You can watch some videos here:
by Prof. Nathan Kutz.
Abhik Saha
Abhik Saha am 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.

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

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

Start Hunting!

Translated by