Why doesn't M2 function get called? Could somebody please explain what this error is and how to fix it?

I have function for M2

function M2 =fRK4M2(M1,M2)
M1=10;
M2=0;
mu_2=10^-3;
K1= 10^-4;
K2=5*10^-4;

M2 = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);

This is my M1 function

function M1 =fRK4M1(M1,M2,M3,O,P)  

  delta=50;
  K1= 10^-4;
  Ko=0.1;
  n=3;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  K2=5*10^-4;
  K3=10^-3;
  gamma=75;
 
  M1=(delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2.*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1)); 
    

And this is my main function

clc;clear;
  
%input for time
t(1)=0;
dt=0.1; %time interval
t=0:dt:100; %time span

%input empty array
T=zeros(length(t),1); %empty array for t
M1=zeros(length(t),1); %empty array for M1
M2=zeros(length(t),1); %empty array for M2

for j = 1:length((t))
    T(j+1)=T(j)+dt;
    M1(j+1)= M1(j)+1./(1+exp(-T(j)));
   
    k1M2 = dt*fRK4M2(M2(j),M1(j));
    k2M2 = dt*fRK4M2(M2(j)+k1M2/2,M1(j)+k1M2/2);
    k3M2 = dt*fRK4M2(M2(j)+k2M2/2,M1(j)+k2M2/2);
    k4M2 = dt*fRK4M2(M2(j)+k3M2,M1(j)+k3M2/2);
    M2(j+1) = M2(j)+1/6*(k1M2+2*k2M2+2*k3M2+k4M2);
    
    k1M1= dt*fRK4M1(M1(j),M2(j));
    k2M1= dt*fRK4M1(M2(j)+k1M2/2,M1(j)+k1M1/2);
    k3M1= dt*fRK4M1(M2(j)+k2M2/2,M1(j)+k2M1/2);
    k4M1= dt*fRK4M1(M2(j)+k3M2/2,M1(j)+k3M1/2);
    M1(j+1) = M1(j)+(1/6*(k1M1+(2*k2M1)+(2*k3M1)+k4M1));
    end

the results I get are not as desired, please someone tell me my fault adn how to fix it. Thanks

[Torsten]

So the task of your code is to solve the system of differential equations

dM1/dt = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2.*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1))
dM2/dt = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2)

?

Since your code has so many errors, I suggest you use ODE45 for the solution.

[cindyawati cindyawati]

thank you for response @Torsten yes my code to solve the system of differential equation. If I use Runge Kutta orde 4th is it possible?

Answer 1

[Lakshya]

Hi,

The code you provided initializes M1 and M2 to specific values (M1 = 10 and M2 = 0) within the fRK4M2 function, which might override the values passed as arguments. If you intend to use the passed arguments instead, you can remove those assignments from the fRK4M2 function.

Hope this helps.

[cindyawati cindyawati]

Hallo @Lakshya thank you for your response, so I must remove M1=10 and M2= 0 in FRK4M2 function? If I remove that , my FRK4M2 function it doesn't run

[Lakshya]

I have saved FRK4M2, FRK4M1 and main functions in seperate files, and now when I'm running the main function, there is the below error in main function. You are passing only 2 arguments in FRK4M1 instead of 5.

you need to provide the missing input arguments when calling the fRK4M1.

M3, O, and P are not defined/initialized.

Not enough input arguments.

Error in fRK4M1 (line 14)
  M1=(delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2.*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1));

Error in answer (line 18)
    k1M1= dt*fRK4M1(M1(j),M2(j));

[Image Analyst]

The error message is fairly explicit -- you're not passing in all the arguments that it needs. You define fRK4M1 as needing 5 input arguments. You are passing it only two. Why? Why are you not passing it 5?

[Lakshya]

Yaa I was telling cindyawati this exact point only ( it seems, it looked like a question rather than a statement ). I am just curious as to why the OP didnt get this error.

[cindyawati cindyawati]

Hello @Lakshya I already input argument for M3,O, and P but my plot for M2, O and P not as expected ( I expected get like sigmoid graph). Why doesnt my graph show a sigmoid graph even though I have include in the main function? This is my main function

clc;clear;
  
%input for time
t(1)=0;
dt=0.1; %time interval
t=0:dt:100; %time span

%input empty array
T=zeros(length(t),1); %empty array for t
M1=zeros(length(t),1); %empty array for M1
M2=zeros(length(t),1); %empty array for M2
M3=zeros(length(t),1); %empty array for M3
O=zeros(length(t),1); %empty array for O
P=zeros(length(t),1); %empty array for P

 %initial condition
    M1 =10;
    M2 = 0;
    M3 = 0;
    O =0;
    P =0;

for j = 1:length((t))
    T(j+1)=T(j)+dt;
    M1(j+1)= M1(j)+1./(1+exp(-T(j)));
    
    k1M2 = dt*fRK4M2(M2(j),M1(j));
    k2M2 = dt*fRK4M2(M2(j)+k1M2/2,M1(j)+k1M2/2);
    k3M2 = dt*fRK4M2(M2(j)+k2M2/2,M1(j)+k2M2/2);
    k4M2 = dt*fRK4M2(M2(j)+k3M2,M1(j)+k3M2/2);
    M2(j+1) = M2(j)+1/6*(k1M2+2*k2M2+2*k3M2+k4M2);
    
    k1M3 = dt*fRK4M3(M1(j),M3(j));
    k2M3 = dt*fRK4M3(M3(j)+k1M3/2,M1(j)+k1M3/2);
    k3M3 = dt*fRK4M3(M3(j)+k2M3/2,M1(j)+k2M3/2);
    k4M3 = dt*fRK4M3(M3(j)+k3M3,M1(j)+k3M3);
    M3(j+1) = M3(j)+1/6*(k1M3+2*k2M3+2*k3M3+k4M3);
    
    k1O = dt*fRK4O(M1(j),O(j));
    k2O = dt*fRK4O(O(j)+k1O/2,M1(j)+k1O/2);
    k3O = dt*fRK4O(O(j)+k2O/2,M1(j)+k2O/2);
    k4O = dt*fRK4O(O(j)+k3O,M1(j)+k3O);
    O(j+1) = O(j)+1/6*(k1O+2*k2O+2*k3O+k4O);
    
    k1P = dt*fRK4P(P(j),M1(j));
    k2P = dt*fRK4P(P(j)+k1P/2,M1(j)+k1P/2);
    k3P = dt*fRK4P(P(j)+k2P/2,M1(j)+k2P/2);
    k4P = dt*fRK4P(P(j)+k3P,M1(j)+k3P/2);
    P(j+1) = P(j)+1/6*(k1P+2*k2P+2*k3P+k4P);
    
    k1M1= dt*fRK4M1(M1(j),M2(j),M3(j),O(j),P(j));
    k2M1= dt*fRK4M1(M2(j)+k1M2/2,M3(j)+k1M3/2,O(j)+k1O/2,P(j)+k1P/2,M1(j)+k1M1/2);
    k3M1= dt*fRK4M1(M2(j)+k2M2/2,M3(j)+k2M3/2,O(j)+k2O/2,P(j)+k2P/2,M1(j)+k2M1/2);
    k4M1= dt*fRK4M1(M2(j)+k3M2/2,M3(j)+k3M3/2,O(j)+k3O/2,P(j)+k3P/2,M1(j)+k3M1/2);
    M1(j+1) = M1(j)+(1/6*(k1M1+(2*k2M1)+(2*k3M1)+k4M1));
    end

Unrecognized function or variable 'fRK4M2'.

  
figure;
plot (T,M1,'r','Linewidth',3)
xlabel('time')
ylabel('M1')

figure 
plot (T,M2,'b','Linewidth',3)
xlabel('time')
ylabel('M2')

figure 
plot (T,M3,'g','Linewidth',3)
xlabel('time')
ylabel('M3')

figure 
plot (T,O,'b','Linewidth',5)
xlabel('time')
ylabel('O')

figure 
plot (T,P,'r','Linewidth',5)
xlabel('time')
ylabel('P')

This is my result plot. Why my O and P (P graph same like O graph) graph no changes

[cindyawati cindyawati]

This is my function for M2,M3, O, and P @Lakshya

function M2 =fRK4M2(M2,M1)
M1=10;
%M2=0;
mu_2=10^-3;
K1= 10^-4;
K2=5*10^-4;     

M2 = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);

M3 function

function M3 =fRK4M3(M3,M1)
%M1=10;
M2=0;
%M3=0;
mu_3=10^-3;
K2=5*10^-4;
K3=10^-3;

M3=(K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);

O function

function O =fRK4O(O,M1)
%M1=10;
M3=0;
O=0;
Ko=0.1;
mu_o=10^-4;
K3=10^-3;


O  = (K3*M1*M3)-(Ko*M1*O)-(mu_o*O);

[Torsten]

The equations you are trying to solve are coupled. Thus you must use Runge-Kutta in a vectorized form - you can't update M1, M2, M3, O and P "one after the other" as you do in your code, but simultaneously. With your abilities to program and your numerical background, you only have a chance to solve the system if you transfer your problem to a ready-to-use integrator. MATLAB has many of them, e.g. ODE45 as already mentionned.

[cindyawati cindyawati]

Hello @Torsten I try using ODE45 and I still confused but my code give error like this

function CM1= mymode (t,M1)

M1= 10;
M2 =0;
delta=50;
K1= 10^-4;
Ko=0.1;
n=3;
Oa=10;
Pa=100;
mu_1=10^-3;
K2=5*10^-4;
K3=10^-3;
gamma=75;
mu_2=10^-3;

CM1= zeros(2,1);
CM1(1)= (delta*M1*(1-(M1./gamma))-2*K1.*M1*M1-M1.*(K2.*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1)); 
CM1(2)=(K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);

end

[t,M1]=ode45(@mymode, [0 100],[0 1]);


Error: File: mymode.m Line: 22 Column: 1
This statement is not inside any function.
 (It follows the END that terminates the definition of the function "mymode".)

[Stephen23]

Function definitions in a script must occur after all other code in the script:

https://www.mathworks.com/help/matlab/matlab_prog/local-functions-in-scripts.html

So for your code:

[t,M1] = ode45(@mymode, [0,100],[0,1])

t = 85×1

1.0e+00 * 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

M1 = 85×2

0 1.0000 0.0001 1.0000 0.0001 1.0000 0.0002 1.0000 0.0002 1.0000 0.0005 1.0000 0.0007 1.0000 0.0010 1.0000 0.0012 1.0000 0.0025 1.0000

function CM1 = mymode (t,M1)
M1= 10; % this overwrite the input value
M2 = 0;
M3 = 1; % you did not define this value
O  = 1; % you did not define this value
delta=50;
K1= 10^-4;
Ko=0.1;
n=3;
Oa=10;
Pa=100;
mu_1=10^-3;
K2=5*10^-4;
K3=10^-3;
gamma=75;
mu_2=10^-3;
CM1= zeros(2,1);
CM1(1) = (delta*M1*(1-(M1./gamma))-2*K1.*M1*M1-M1.*(K2.*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1));
CM1(2) = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);
end

[cindyawati cindyawati]

Why my graph cannot show like sigmoid graph whereas I already input the sigmoid function into the code?

[Torsten]

If you mean the code with ODE45, it's most probably because you solve the wrong differential equation. But we cannot tell.

Where do you input a sigmoid function ?

[cindyawati cindyawati]

No, I mean the code with Runge-Kutta orde 4. I input sigmoid function inside for looping before runge kutta. This is my code

clc;clear;
  
%input for time
t(1)=0;
dt=0.1; %time interval
t=0:dt:100; %time span

%input empty array
T=zeros(length(t),1); %empty array for t
M1=zeros(length(t),1); %empty array for M1
M2=zeros(length(t),1); %empty array for M2

for j = 1:length((t))
    T(j+1)=T(j)+dt;
    M1(j+1)= M1(j)+1./(1+exp(-T(j))); %this is sigmoid function
   
    k1M2 = dt*fRK4M2(M2(j),M1(j));
    k2M2 = dt*fRK4M2(M2(j)+k1M2/2,M1(j)+k1M2/2);
    k3M2 = dt*fRK4M2(M2(j)+k2M2/2,M1(j)+k2M2/2);
    k4M2 = dt*fRK4M2(M2(j)+k3M2,M1(j)+k3M2/2);
    M2(j+1) = M2(j)+1/6*(k1M2+2*k2M2+2*k3M2+k4M2);
    
    k1M1= dt*fRK4M1(M1(j),M2(j));
    k2M1= dt*fRK4M1(M2(j)+k1M2/2,M1(j)+k1M1/2);
    k3M1= dt*fRK4M1(M2(j)+k2M2/2,M1(j)+k2M1/2);
    k4M1= dt*fRK4M1(M2(j)+k3M2/2,M1(j)+k3M1/2);
    M1(j+1) = M1(j)+(1/6*(k1M1+(2*k2M1)+(2*k3M1)+k4M1));
    end

Unrecognized function or variable 'fRK4M2'.

This my function for M1, M2, M3, O, and P

function M2 =fRK4M2(M2,M1)
M1=10;
%M2=0;
mu_2=10^-3;
K1= 10^-4;
K2=5*10^-4;     

M2 = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);

function M3 =fRK4M3(M3,M1)
%M1=10;
M2=0;
%M3=0;
mu_3=10^-3;
K2=5*10^-4;
K3=10^-3;

M3=(K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);


function O =fRK4O(O,M1)
%M1=10;
M3=0;
O=0;
Ko=0.1;
mu_o=10^-4;
K3=10^-3;


O  = (K3*M1*M3)-(Ko*M1*O)-(mu_o*O);

function P =fRK4P(P,M1)
%M1=10;
O=0;
%P=0;
Ko=0.1;
mu_p= 10^-5;

P  = (Ko*M1*O)-(mu_p*P);
 
function M1 =fRK4M1(M1,M2,M3,O,P)  

  delta=50;
  K1= 10^-4;
  Ko=0.1;
  n=3;
  Oa=10;
  Pa=100;
  mu_1=10^-3;
  K2=5*10^-4;
  K3=10^-3;
  gamma=75;
 
  M1=(delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2.*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1)); 
    

[Torsten]

I already told you that your Runge-Kutta code is incorrect in many respects.

Concerning the sigmoid function:

You update M1 as

M1(j+1)= M1(j)+1./(1+exp(-T(j))); %this is sigmoid function

and then overwrite this setting at the end of the loop by

M1(j+1) = M1(j)+(1/6*(k1M1+(2*k2M1)+(2*k3M1)+k4M1));

What is the motivation for this ?

And there are many, many other examples of this kind in the Runge-Kutta code.

So I can only repeat: Set up your problem using an ODE solver. This way, you have at least the guarantee that the numerical method integrates your ODE system numerically correct. Whether you set up the system correctly is a different question.

[cindyawati cindyawati]

Ok, I know where I went wrong when I used the 4th order Runge kutta formula. I tried removing the sigmoid function and my program is still running. But the Y-axis for M3 has a negative value and the O and P graphs remain flat. Why?

[cindyawati cindyawati]

Hello @Torsten I wanna ask If I use ODE45 for my problem. I must make 5 function (M1,M2, M3,O,and P) for ODE45 in separate file or I make one file include 5 function? Btw I already try solve my code using ODE45 but i get error like this. How to fix it? Thank you

function CM1 = mymode (t,M1)
M1= 10; 
M2 = 0;
M3 = 0; 
O  = 0; 
P=0;

delta=50;
K1= 10^-4;
Ko=0.1;
n=3;
Oa=10;
Pa=100;
mu_1=10^-3;
K2=5*10^-4;
K3=10^-3;
gamma=75;
mu_2=10^-3;
mu_3=10^-3;
mu_o=10^-4;
mu_p= 10^-5;

CM1= zeros(2,1);
CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1));
CM1(2) = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);
CM1(3) = (K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);
CM1(4) = (K3*M1*M3)-(Ko*M1*O)-(mu_o*O);
CM1(5) = (Ko*M1*O)-(mu_p*P);

[t,M1] = ode45(mymode, [0,100],[0,1])
plot (t,M1)
end


>> mymode
Out of memory. The likely cause is an infinite recursion within the program.

Error in mymode (line 31)
[t,M1] = ode45(mymode, [0,100],[0,1])