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

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]);

[Dyuman Joshi]

"Please help some one"

What do you need help with? You should specify that in your question.

[mks]

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]

Sorry typing mistake badin=basin

Answer 1

[Sam Chak]

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

[mks]

here x,y,z are population density. you have chosen y as negative but it is not possible. For the parameters set we are getting three coexistence equilibrium point. You should change the initial condition.

[Sam Chak]

Have you tested your supplied initial condition in the code?

I didn't know about the nature of your ODEs. The stable equilibrium point was found by accident.

[Sam Chak]

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]

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

[Sam Chak]

@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.