Bifurcation diagram for a system of 3 non-linear ODE

2 次查看（过去 30 天）
Giovanni Viciconte
评论： Delfim Joao 2019-11-29
Hi everybody. I have a little problem to plot a bifurcation diagram for a system of 3 ODE. This is my system. The parameter is "da"
function dAdt = SistD( t,A )
global da mua mub nua nub nua_s etab xia xib
dAdt = zeros(3,1);
e=1e-2;
la=2.01+71*da;
la=la*e;
lb=-1.85+76.4i+(66.7+9.35i)*da;
lb=lb*e;
dAdt(1)=la*A(1) -mua*A(1)*abs(A(1))^2 -nua*A(1)*abs(A(2))^2 -nua_s*A(1)*abs(A(3))^2 -xia*conj(A(1))*conj(A(3))*A(2);
dAdt(2)=lb*A(2) -mub*A(2)*abs(A(2))^2 -nub*A(2)*abs(A(3))^2 -etab*A(2)*abs(A(1))^2 -xib*A(3)*(A(1))^2;
dAdt(3)=lb*A(3) -mub*A(3)*abs(A(3))^2 -nub*A(3)*abs(A(2))^2 -etab*A(3)*abs(A(1))^2 -xib*A(2)*(conj(A(1)))^2;
end
In order to obtain the bifurcation diagram I solve the system with ode45. In this case I use the sum of the modulus of the three variables "phi" as "global variable" to understand the system behavior. I use the "stationary" value of the solution(the value that the solution reach for t>>0),the last 19 values of the solution, to plot the diagram.
clear variables
close all
clc
global da mua mub nua nub nua_s etab xia xib
e=1e-2;
%%Coeff A
mua=3.11;
mua=mua*e;
nua=6.88-1.11i;
nua=nua*e;
nua_s=6.88+1.11i;
nua_s=nua_s*e;
xia=4.57;
xia=xia*e;
%%Coeff B
mub=2.42+0.0321i;
mub=mub*e;
nub=3.13-0.816i;
nub=nub*e;
etab=0.955-3.47i;
etab=etab*e;
xib=1.62-1.36i;
xib=xib*e;
%%Calc
Rec=121.1;
eps=1e-1;
TF=500;
tspan=[0 700];
PHI=zeros(32,20);
options = odeset('RelTol',1E-12,'AbsTol',1E-15,'MaxStep',0.01);
%%Bif
for Re=115:147
A0=[3 3 3]';
da=(1/Rec-1/Re)*1/eps^2;
[t,L]=ode45(@SistD,tspan,A0,options);
TF=max(size(L));
phi=abs(L(:,1))+abs(L(:,2))+abs(L(:,3));
PHI(Re-114,:)=phi(TF-19:TF)';
end
plot(115:1:147,PHI)
But I do not have the good solution. My questions are 1) Is my method right for this kind of non-linear system ? 2) The fact that i use conj function does affect the effectiveness of ode45? Thank you for your time and help

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Nonlinear Analysis 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by