Why my y value went so high, whereas I using small value?
조회 수: 2 (최근 30일)
이전 댓글 표시
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1411409/image.png)
댓글 수: 2
Image Analyst
2023년 6월 15일
People may wait to answer until you attach the fRK4M2 variable or function so that they can run your code.
채택된 답변
Steven Lord
2023년 6월 15일
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
cindyawati cindyawati
2023년 6월 15일
편집: cindyawati cindyawati
2023년 6월 15일
추가 답변 (0개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!