Why my y value went so high, whereas I using small value?
1 view (last 30 days)
Show older comments
cindyawati cindyawati
on 15 Jun 2023
Edited: cindyawati cindyawati
on 15 Jun 2023
I am plotting M1 with RK4 function. I am just using small value of constant but the y value went so high
%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
for j = 1:length((t))
T(j+1)=T(j)+dt;
M1(j+1)= M1(j)+1./(1+exp(-T(j)));
k1M1= dt*fRK4M1(M2(j),M3(j),O(j),P(j),M1(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')
This is my RK function for M1
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;
mu_2=10^-3;
mu_3=10^-3;
mu_o=10^-4;
mu_p= 10^-5;
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 plot
2 Comments
Image Analyst
on 15 Jun 2023
People may wait to answer until you attach the fRK4M2 variable or function so that they can run your code.
Accepted Answer
Steven Lord
on 15 Jun 2023
Look at the way you're calling your function in the loop:
k1M1= dt*fRK4M1(M2(j),M3(j),O(j),P(j),M1(j));
and the way you've defined it:
function M1 =fRK4M1(M1,M2,M3,O,P)
Did you mean for M1 to be the first input argument in the definition but the last in the call? That smells like a bug.
There's also nothing in your code that changes any element of M2, M3, O, or P from their initial values of 0 so you could save time by simply omitting them from the calls (replacing them with 0.)
1 Comment
cindyawati cindyawati
on 15 Jun 2023
Edited: cindyawati cindyawati
on 15 Jun 2023
More Answers (0)
See Also
Categories
Find more on 2-D and 3-D Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!