MATLAB: Ode solver results not as expected

ode solver wrong results

Hello,

I am trying to solve differential equations numerically, but the results I obtain are not as expected. I should be able to plot an increasing exponential as a solution, but instead I get only a zero vector.

This is my differential equation:

dY/dt=(ßY^n)/(K^n+Y)^n)-αY

This is the script I am using:

%Plot of Simple Gene Regulation and NAR

clear all;

tspan = [0 100]; Y0 = 0; X0 = 0; beta = 8; alfa = 0.101; n = 4; k = 9; N = 1;

Yss = beta/alfa; %Xss = k*((beta/alfa)*k)^(1/(n+1)); Xss = k;

[t1,Y] = ode23(@(t,Y) (beta/(1+(k/Y)^N)-alfa*Y), tspan, Y0); %Simple Reg [t2,X] = ode45(@(t,X) (beta/(1+(X/k)^n)-alfa*X), tspan, X0); %NAR

%plot(t1,Y./Yss,t2,X./Xss); plot(t1,Y); grid on; title('Comparison between Simple Gene Regulation and NAR'); xlabel('Time t'); ylabel('Normalized Steady States'); %legend('Simple Regulation','NAR'); hold on;

The answers I get are good if I use simulink to calculate the solutions only:

I should be able to obtain the same plots by using the command line ode solvers, but there is something wrong I might be missing. Any help would be greatly appreciatted. Thank you very much.

Best regards,

Alejandro

Best Answer

It looks like expected:

dY/dt = (ß * Y^n) / (K^n + Y^n) - α * Y

If the initial value of Y is 0, the derivative is 0 also:

dY/dt = (ß * 0^n) / (K^n + 0^n) - α * 0

In consequence the trajectory remains at the 0 level.

The Y^n in the numerator is useful to avoid a division by 0.

Related Solutions

MATLAB: 4th order Runge-Kutta code that can solve for several intial conditions

Here's a version:

% Initial conditions
t0=0;
x0=0:0.2:1;
y0=1:0.2:2;
dt=0.1;
tz=1;
t_range= t0:dt:tz;
X = zeros(numel(x0),numel(t_range));
Y = zeros(numel(y0),numel(t_range));
x_rk=zeros(1,numel(t_range));
y_rk= zeros(1,numel(t_range));
for n = 1:numel(x0)
    x_rk(1)=x0(n);
    y_rk(1)=y0(n);
    for i=1:length(t_range)-1
        [x_rk(i+1), y_rk(i+1)]=rk_method(x_rk(i),y_rk(i),dt);
    end
    X(n,:) = x_rk;
    Y(n,:) = y_rk;
    %% PLots
    subplot(numel(x0)/2,2,n)
    plot (t_range, x_rk,'b-o',t_range, y_rk,'y-o','MarkerSize',5);
    grid
    title([num2str(x0(n)),'   ',num2str(y0(n))])
    legend('x','y')
end
disp('X = '), disp(X)
disp('Y = '), disp(Y)
%% Functions
function [dxdt, dydt]= f(x,y)
dxdt =  5*x-3*y;
dydt = -6*x+2*y;
end 
function [x1, y1]=rk_method(x0,y0,dt)
[k11, k12]=f(x0,y0);     
[k21, k22]=f(x0+0.5*dt*k11,y0+0.5*dt*k12); 
[k31, k32]=f(x0+0.5*dt*k21,y0+0.5*dt*k22); 
[k41, k42]=f(x0+dt*k31, y0+dt*k32); 
x1=x0+dt*((k11/6)+((k21+k31)/3)+(k41/6));
y1=y0+dt*((k12/6)+((k22+k32)/3)+(k42/6));
end

MATLAB: Plotting two variables from a loop

% Constants
g=9.81;
m=1.5;
A=0.028; 
Cd=0.8;
rho=1.225;
B=(rho*Cd*A)/2;
Dt=1e-6;
% Initial values
x=0;
y=6.4008; %21 ft
V=12.5;
alpha=90;
Vx=V*sind(alpha);
Vy=V*cosd(alpha);
t=0;
X = zeros([],1) ; 
Y = zeros([],1) ; 
for k=1:1e9
    
    if y<0
        break
    end
    
    ax= -(B/m)*(V*Vx);
    ay=(B/m)*(V*Vy)-g;
    x=x+Vx*Dt+0.5*ax*(Dt^2);
    y=y+Vy*Dt+0.5*ay*(Dt^2);
    Vx=Vx+ax*Dt;
    Vy=Vy+ay*Dt;
    V=sqrt((Vx^2)+(Vy^2));
    
    t=t+Dt;
    X(k) = x ;
    Y(k) = y ; 
end
plot(X,Y)

Related Question

Solving a linear system of 1st order ODEs using the Forward Euler method and ode45