I want to isocline induced step function

x'=rx(1-x/k)-axy/(1+quM)

y'=(eaxy)/(1+quM)-my

Where u=0 if x<=pi_0; u=sx/(sx+F)+w(sx-F)/(b(sx+F)M), if pi_0 <x<pi_1; u=1 if x>=pi_1

where pi_0=(wF)/(s(w+bM))

pi_1=F(w+bM)/(sw)

Answer 1

[Mrutyunjaya Hiremath]

• To simulate the isocline-induced step function differential equations, we can use MATLAB's ODE solver `ode45`.
• First define the system of differential equations as a function and then solve it numerically.

• 1. Define the system of differential equations:

function dydt = isocline_system(t, y, r, k, a, e, m, F, s, w, b, M, sx)
    x = y(1);
    u = isocline_input(x, F, s, w, b, M, sx);

    dxdt = r * x * (1 - x / k) - a * x * y / (1 + q(u) * u * M);
    dydt = (e * a * x * y) / (1 + q(u) * u * M) - m * y;
    dydt = [dxdt; dydt];
end

function u = isocline_input(x, F, s, w, b, M, sx)
    pi0 = (w * F) / (s * (w + b * M));
    pi1 = F * (w + b * M) / (s * w);

    if x <= pi0
        u = 0;
    elseif x >= pi1
        u = 1;
    else
        u = sx / (sx + F) + w * (sx - F) / (b * (sx + F) * M);
    end
end

function q_u = q(u)
    q_u = 1 / (1 - u);
end

• 2. Solve the differential equations numerically using `ode45`:

% Parameters
r = ...; % Define the value of r
k = ...; % Define the value of k
a = ...; % Define the value of a
e = ...; % Define the value of e
m = ...; % Define the value of m
F = ...; % Define the value of F
s = ...; % Define the value of s
w = ...; % Define the value of w
b = ...; % Define the value of b
M = ...; % Define the value of M
sx = ...; % Define the value of sx

% Initial conditions
x0 = ...; % Define the initial value of x
y0 = ...; % Define the initial value of y
initial_conditions = [x0; y0];

% Time span for simulation
tspan = ...; % Define the time span (e.g., [0, 10])

% Solve the system of differential equations
[t, y] = ode45(@(t, y) isocline_system(t, y, r, k, a, e, m, F, s, w, b, M, sx), tspan, initial_conditions);

% Extract the solution
x_sol = y(:, 1);
y_sol = y(:, 2);

• Replace the `...` with the specific values of the parameters and initial conditions you want to use in your simulation. The `ode45` solver will numerically integrate the system of differential equations over the given time span `tspan`, and the solutions for `x` and `y` will be stored in `x_sol` and `y_sol`, respectively.

[mks]

% Parameters

r =0.1; % Define the value of r

k = 50; % Define the value of k

a =0.01; % Define the value of a

e = 0.5; % Define the value of e

m = 0.05; % Define the value of m

F = 1; % Define the value of F

%s =0.1; % Define the value of s

w = 0.1; % Define the value of w

b = 0.001; % Define the value of b

M = 50; % Define the value of M

s =0.1; % Define the value of sx

% Initial conditions

x0 =2; % Define the initial value of x

y0 =1; % Define the initial value of y

initial_conditions = [x0; y0];

% Time span for simulation

tspan = [0, 100]; % Define the time span (e.g., [0, 10])

% Solve the system of differential equations

[t, y] = ode45(@(t, y) kkk(t, y, r, k, a, e, m, F, s, w, b, M, x), tspan, initial_conditions);

% Extract the solution

x_sol = y(:, 1);

y_sol = y(:, 2);

function dydt = kkk(t,x,y, r, k, a, e, m, F, s, w, b, M)

x = y(1);

u = isocline_input(x, F, s, w, b, M);

dxdt = r * x * (1 - x / k) - a * x * y / (1 + q * u * M);

dydt = (e * a * x * y) / (1 + q * u * M) - m * y;

dydt = [dxdt; dydt];

end

function u = isocline_input(x, F,w, b, M, s)

pi0 = (w * F) / (s * (w + b * M));

pi1 = F * (w + b * M) / (s * w);

if x <= pi0

u = 0;

elseif x >= pi1

u = 1;

else

u = (s*x) / (s*x + F) + w * (s*x - F) / (b * (s*x + F) * M);

end

end

I want to plot (x,y) where is my fault please help

Answer 2

[Sam Chak]

Hi @MD KAUSAR SK

The above code seems to be generated by ChatGPT and it has flaws. You shouldn't directly use or blindly trust someone who simply gave you an untested code without showing how it work.

ic     = [2; 1];    % initial values
tspan  = [0 100];   % sim time
[t, z] = ode45(@odefcn, tspan, ic);

x      = z(:, 1);
y      = z(:, 2);
subplot(2,1,1)
plot(t, x), grid on, title('x(t)'), xlabel('t')
subplot(2,1,2)
plot(t, y), grid on, title('y(t)'), xlabel('t')

function dzdt = odefcn(t, z)
    % parameters
    r = 0.1; 
    k = 50;  
    a = 0.01;
    q = 0.03;   % <--- change this value
    M = 50;
    e = 0.5;
    m = 0.05;
    w = 0.1;
    F = 1.613;  % <--- change this value
    s = 0.1;
    b = 0.001;
        
    x = z(1);
    y = z(2);
    
    % system input
    pi0 = (w*F)/(s*(w + b*M));
    pi1 = F*(w + b*M)/(s*w);
    
    if x <= pi0
        u = 0;
    elseif x >= pi1
        u = 1;
    else
        u = (s*x)/(s*x + F) + w*(s*x - F)/(b*(s*x + F)*M);
    end
    
    % dynamics
    dxdt = r*x*(1 - x/k) - (a*x*y)/(1 + q*u*M);
    dydt = (e*a*x*y)/(1 + q*u*M) - m*y;
    dzdt = [dxdt; 
            dydt];
end