Why my plot cannot show sigmoid graph? Whereas I already input to the function

Question

Hello I am confused with my code in MATLAB. I am use Runge-Kutta orde 4th. Why my plot cannot show sigmoid whereas I already input the sigmoid function into my code. Please somebody can explain how to fix it. 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
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))); % sigmoid eq
    
    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(O(j),M1(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
  
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 function for M1, M2, M3, O, and P in separate file

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

my result like this my plot for M1,M2, and M3 and for O and P like this my plot for O and P

I can only guess the math is wrong. What happens is that `O` is always zero. Your function `fRK4O` will always return zero if `O` is zero, because it returns `O = (K3*M1*M3)-(Ko*M1*O)-(mu_o*O);`. The first term is always zero (`M3=0`) and the other two are zero if `O=0`. But this is the function you use to update `O`, there fore it will always be zero. — Ander Biguri, Jun 19 '23 at 13:10
Sorry, I mean the value of M2, M3, O, and P is always 0 so I cant replace value of M2,M3, O and P — Cindy, Jun 19 '23 at 13:36
I am not suggesting you replace it. Perhaps I misunderstood your question, I thought you were asking why `O` is zero in your plot. I explained why. Can you please explain what is really the question you have? What of what you are showing here is not doing what you want, and what do you expect it to do instead? — Ander Biguri, Jun 19 '23 at 13:37
I want to solve a differential equation using a 4th order Runge Kutta. I have written the function I used above. I want to get a graph which is sigmoid in shape. The function I wrote above when I run it doesn't give the result I want. Even though I've included the sigmoid function and hope my graph looks like a sigmoid — Cindy, Jun 19 '23 at 13:49
Which of the plots should look like a sigmod? all of them? which one you say its a sigmoid? I think in general your equations are wrong, so its not a programming question. But you have not provided where your equations are frmo/which differential equation you are trying to solve, so we can't help. For example, I explained why in your code, `O` will always be `O=0`, so I suspect your math to compute `O` is wrong. But I don't know what `O` is supposed to be, you have not explained. — Ander Biguri, Jun 19 '23 at 13:52
Look at the line `k2M2 = dt*fRK4M2(M2(j)+k1M2/2,M1(j)+k1M2/2);` and contemplate what could be wrong about it, just from the variable names and their intended meaning. The stage computations after it have the same problem, this is just the first instance. — Lutz Lehmann, Jun 19 '23 at 16:06
it should be k2M2=dt*fRK4M2(M1(j)+k1M2/2+M2(j)+k1M2/2)? @LutzLehmann — Cindy, Jun 20 '23 at 07:09
No, this is the same error only in another place. The updates for M1 are the k1M1, k2M1, k3M1 while the updates for M2 are the k1M2, k2M2, k3M2 etc. Updates for O are k1O, k2O,... Of course these need all to be computed first in the sequence of computation before they can be inserted. In explicit methods this is well possible in an organized manner. /// It would help to understand your code better and for a longer time if the outputs had different names than the inputs. So the slope function for O would have a result of dO or dotO or Odot or whatever you find intuitive to mark a derivative. — Lutz Lehmann, Jun 20 '23 at 07:30
@AnderBiguri Yes, I want all of the plot look like a sigmoid. This is the formula I wanted to finish with a 4th order Runge Kutta.(https://drive.google.com/file/d/1ab8IJ5phW1C9nHiyVAvTDen-Oqzl0urn/view?usp=sharing) but my function is couple each other — Cindy, Jun 20 '23 at 07:41
@LutzLehmann do you mean I have to create a new variable to store the calculation results? Like this M2new= M2(j+1)? — Cindy, Jun 20 '23 at 08:05
With the last name comment I mean that you should declare your functions as `function dotO =fRK4O(O,M1)`. As it is not specific to the RK4 method, you could also name it as `function dotO = derivs_O(O,M1)`, or `O_prime` or so. // You are computing `M1(j+1)` twice, what is the real version? You can not compute `k2M3` if you do not have `k1M1`. // You should check for argument order, like in `k1M3 = dt*fRK4M3(M1(j),M3(j));` or `k1M1= dt*fRK4M1(M1(j),M2(j),M3(j),O(j),P(j));` which are different to the lines following them. — Lutz Lehmann, Jun 20 '23 at 08:18
For the non-trivial error of not preserving the coupled nature of the system see https://stackoverflow.com/questions/46654283/runge-kutta-for-coupled-odes/53996811#53996811, https://stackoverflow.com/questions/58255653/lotka-volterra-with-runge-kutta-not-desired-precision/58258912#58258912, https://stackoverflow.com/questions/64867616/rk4-integrator-of-a-nonlinear-second-order-ode-in-python/64873723#64873723 — Lutz Lehmann, Jun 20 '23 at 08:24
sorry @LutzLehmann for the function I already declare in separate files like this 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); but that didn't work. // Which order is correct ? k1M1= dt*fRK4M3 (M1(j), M3(j)); or k1M1= dt*fRK4M3 (M3(j), M1(j)) // In which part shoul I check whether in the function section (like function M1) or in the main function? thnk — Cindy, Jun 20 '23 at 08:31
You are the programmer, you decide what the correct argument order is. // The names are mostly not important for the execution of the program, but for avoiding programming errors and to ease debugging and maintaining. // If the math is not reflected correctly in the program, then the results will not reflect the math. — Lutz Lehmann, Jun 20 '23 at 09:09

Why my plot cannot show sigmoid graph? Whereas I already input to the function

0 Answers0