I want to plot Basin of attraction. I have written a few line of code given below.Please help some one

Question

mks 2023-8-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2004982-i-want-to-plot-basin-of-attraction-i-have-written-a-few-line-of-code-given-below-please-help-some

评论： Sam Chak 2023-8-6

clear all
clc;
r1=0.8;  K1=8;c1=0.12; alpha1=0.505;alpha2=0.5;gamma=0.25;d1=0.2;
r2=0.5;c2=0.22; K2=8;k=1;beta=0.27;lambda=1.07;d=0.3;h=1.2;
f=r1*x(1)-(c1*(x(1))^2)/(k1+alpha1*x(2));
Unrecognized function or variable 'x'.
g=(r2*(x(2))^2)/(1+k*x(3))-(c2*(x(2))^2)/(k2+alpha2*x(1)*x(2))-(beta*(x(2))^2*x(3))/(beta*h*(x(2))^2+x(2)+gamma);...
    h=(lamda*beta*(x(2))^2*x(3))/(beta*h*(x(2))^2+x(2)+gamma)-d*x(3)-d1*x(3)^2;
[x0 y0 z0] = meshgrid(0:0.01:310, 0:0.02:100, 0:0.02:30)
distance = sqrt((x - 306.33859).^2 + (y - 75.150076).^2 + (z -3.29535).^2);
error=0.001;
for i = 1:numel(x0)
    x = x0(i);
    y = y0(i);
    z = z0(i);
    
    if distance <error
        basin(i)= 1;
    end
end 
[t, sol] = ode45(@(t, y) [f; g; h], [0, 1000], [x; y; z]);

3 个评论
显示 1更早的评论隐藏 1更早的评论

mks 2023-8-5

移动：Torsten 2023-8-5

Sir, I can not complete the code . Actually I am trying to plot badin of attraction for my mathematical model. I have written some part of the code . please complete my incomplete code. One equilibrium point is (306.33859,75.150076,3.29535) and parameters value are given the above .for this parameter set the system prpduce three interior equilibrium points out of this one is (306.33859,75.150076,3.29535)

mks 2023-8-5

移动：Torsten 2023-8-5

Sorry typing mistake badin=basin

请先登录，再进行评论。

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

Answer 1

Sam Chak 2023-8-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2004982-i-want-to-plot-basin-of-attraction-i-have-written-a-few-line-of-code-given-below-please-help-some#answer_1283687

编辑：Sam Chak 2023-8-5

在 MATLAB Online 中打开

Hi @MD KAUSAR SK

I'm unsure if the basin of attraction exists or not because a random sampling of the initial conditions in the ranges

,

shows that the trajectories diverge to infinity as time t goes to some finite value. Could you perform a stability analysis on the system? The origin

seems to be the equilibrium.

Update: Found a stable equilibrium near the center of this region. You should be able to find the basin of attraction.

opts = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);

% initial condition

numP = 9;

x0 = linspace( 18.3, 28.3, numP);

y0 = linspace(-13.9, -3.9, numP);

z0 = linspace( 0.2, 10.2, numP);

[cx, cy, cz] = ndgrid(x0, y0, z0);

combs = [cx(:), cy(:), cz(:)];

for j = 1:length(combs)

[t, sol] = ode45(@odefcn, [0, 1056], combs(j, :), opts);

x = sol(:,1);

y = sol(:,2);

z = sol(:,3);

plot3(x, y, z, 'color', '#f99954'), hold on

end

grid on, xlabel('x'), ylabel('y'), zlabel('z')

axis square

axis([18 28 -14 -4 0 10])

% the system

function dxdt = odefcn(t, x)

r1 = 0.8;

k1 = 8;

c1 = 0.12;

alpha1 = 0.505;

alpha2 = 0.5;

gamma = 0.25;

d1 = 0.2;

r2 = 0.5;

c2 = 0.22;

k2 = 8;

k = 1;

beta = 0.27;

lambda = 1.07;

d = 0.3;

h = 1.2;

f = r1*x(1) - (c1*x(1)^2)/(k1 + alpha1*x(2));

g = (r2*x(2)^2)/(1 + k*x(3)) - (c2*x(2)^2)/(k2 + alpha2*x(1)*x(2)) - (beta*(x(2)^2)*x(3))/(beta*h*x(2)^2 + x(2) + gamma);

h = (lambda*beta*(x(2)^2)*x(3))/(beta*h*(x(2)^2) + x(2) + gamma) - d*x(3) - d1*x(3)^2;

dxdt = [f;

g;

h];

end

5 个评论
显示 3更早的评论隐藏 3更早的评论

Sam Chak 2023-8-6

编辑：Sam Chak 2023-8-6

在 MATLAB Online 中打开

Hi @MD KAUSAR SK

I tried using the initial conditions

but the trajectory couldn't even go past

s. Please check if the math model in odefcn(t, x) is the same as your hand-derived model, where you found the equilibrium

. Human errors can occur because the equations are quite long.

Update: The state equation g is updated according to this comment. Still, the trajectory couldn't go past

s.

opts = odeset('RelTol', 1e-4, 'AbsTol', 1e-6);

% initial condition

ic = [306.33859, 75.150076, 3.29535];

[t, sol] = ode45(@odefcn, [0, 0.2], ic, opts);

Warning: Failure at t=1.153514e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-16) at time t.

x = sol(:,1);

y = sol(:,2);

z = sol(:,3);

plot3(x, y, z, 'color', '#f99954'), hold on

grid on, xlabel('x'), ylabel('y'), zlabel('z')

axis square

% the system

function dxdt = odefcn(t, x)

r1 = 0.8;

k1 = 8;

c1 = 0.12;

alpha1 = 0.505;

alpha2 = 0.5;

gamma = 0.25;

d1 = 0.2;

r2 = 0.5;

c2 = 0.22;

k2 = 8;

k = 1;

beta = 0.27;

lambda = 1.07;

d = 0.3;

h = 1.2;

f = r1*x(1) - (c1*x(1)^2)/(k1 + alpha1*x(2));

% g = (r2*x(2)^2)/(1 + k*x(3)) - (c2*x(2)^2)/(k2 + alpha2*x(1)*x(2)) - (beta*(x(2)^2)*x(3))/(beta*h*x(2)^2 + x(2) + gamma);

g = (r2*x(2)^2)/(1 + k*x(3)) - (c2*x(2)^2)/(k2 + alpha2*x(1)) - (beta*(x(2)^2)*x(3))/(beta*h*x(2)^2 + x(2) + gamma);

h = (lambda*beta*(x(2)^2)*x(3))/(beta*h*(x(2)^2) + x(2) + gamma) - d*x(3) - d1*x(3)^2;

dxdt = [f;

g;

h];

end

mks 2023-8-6

Please replace the term k2+alpha2*x(1) instead of k2+alpha2*x(1)*x(2) In the function 'g'

Sam Chak 2023-8-6

@MD KAUSAR SK, Updated g (see above).

Can you click this

button and insert your MATLAB code here? So that it can be tested when the

Run button is clicked.

请先登录，再进行评论。

I want to plot Basin of attraction. I have written a few line of code given below.Please help some one

3 个评论
显示 1更早的评论隐藏 1更早的评论

回答（1 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

Community Treasure Hunt

I want to plot Basin of attraction. I have written a few line of code given below.Please help some one

3 个评论 显示 1更早的评论隐藏 1更早的评论

回答（1 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

Community Treasure Hunt

WeChat

3 个评论
显示 1更早的评论隐藏 1更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论