How to declare j in runge kutta orde 4 function? Why my output doesnt exist?

4 次查看（过去 30 天）

cindyawati cindyawati 2023-5-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1968744-how-to-declare-j-in-runge-kutta-orde-4-function-why-my-output-doesnt-exist

编辑： Torsten 2023-5-20

clc;clear;
function RK4
%input konstanta
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^-4;
gamma=75;
%input initial condition
M1(1)=10;
M2(1)=0;
M3(1)=0;
%input for time
t(1)=0;
dt=0.01; %time interval
t=0:dt:100; %time span
%input empty array
T=zeros(length(t)+1,1); %empty array for t
M1=zeros(length(t)+1,1); %empty array for M1
M2=zeros(length(t)+1,1); %empty array for M2
M3=zeros(length(t)+1,1); %empty array for M3
O=zeros(length(t)+1,1);
P=zeros(length(t)+1,1);
fM1=@(M1,T) M1(j+1)+(dt*(delta*M1(j+1)*(1-(M1(j+1)/gamma))-2*K1*M1(j+1)*M1(j+1)-M1(j+1)*(K2.*M2(j+1))-((Oa-n)*K3*M1(j+1)*M3(j+1))-((Pa-Oa)*Ko*M1(j+1)*O(j+1))-(mu_1*M1(j+1))));
for j = 1:length((t)-1) 
    
    %T(j+1)=T(j)+dt;
    M1 =M1(j)+1./(1+exp(-T(j)));
    
    k1M1 = dt*fM1(M1(j),dt);
    k2M1 = dt*fM1(M1(j)+k1M1/2, dt+k1M1/2);
    k3M1 = dt*fM1(M1(j)+k2M1/2, dt+k2M1/2);
    K4M1 = dt*fM1(M1(j)+k3M1,dt+k3M1);
    M1(j+1) = M1(j)+1/6*(k1M1+2*k2M1+2*k3M1+k4M1);
    T(j+1)=T(j)+dt;
    
   
end
figure 
plot (T,M1,'r','Linewidth',3)
xlabel('time')
ylabel('M1')
end

11 个评论
显示 9更早的评论隐藏 9更早的评论

cindyawati cindyawati 2023-5-20

oke I will try again, sorry if this question bothers you

Torsten 2023-5-20

编辑：Torsten 2023-5-20

If the question would bother me, I would not answer to it.

But I'm always surprised why people like you who have very little knowledge about MATLAB don't use MATLAB's onramp training course first before trying to solve rather complicated programming tasks:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

请先登录，再进行评论。

请先登录，再回答此问题。