I want to solve ODE for a complex systems in a rk4 method

4 次查看（过去 30 天）
mks
mks 2024-3-8
评论： mks 2024-3-9
I know to solve for simple system But i want to solve complex systems using matrix , array in rk4 method. I have a simple rk4 code , now i want to extend to a complex system. You can assume variables are vector .
function [x,y] = rk4th(dydx,xo,xf,yo,h)
x = xo:h:xf ;
y = zeros(1,length(x));
y(1)= yo ;
for i = 1:(length(x)-1)
k_1 = dydx(x(i),y(i));
k_2 = dydx(x(i)+0.5*h,y(i)+0.5*h*k_1);
k_3 = dydx((x(i)+0.5*h),(y(i)+0.5*h*k_2));
k_4 = dydx((x(i)+h),(y(i)+k_3*h));
y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end
dydx = @(x,y) 3.*exp(-x)-0.4*y;
%[x,y] = rk4th(dydx,0,100,-0.5,0.5);
%plot(x,y,'o-');
end

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

采纳的回答

James Tursa
James Tursa 2024-3-9
Just turn all the y( ) stuff into row or column vectors. E.g., suppose we use column vectors, then something like this ...
y = zeros(numel(yo),length(x));
y(:,1) = yo;
for i = 1:(length(x)-1)
k_1 = dydx( x(i) , y(:,i) );
k_2 = dydx( x(i)+0.5*h, y(:,i)+0.5*h*k_1 );
k_3 = dydx( x(i)+0.5*h, y(:,i)+0.5*h*k_2 );
k_4 = dydx( x(i)+ h, y(:,i)+ h*k_3 );
y(:,i+1) = y(:,i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end
  1 个评论
mks
mks 2024-3-9
how can I write the two different system of equation below here;
mu=0.0001;h=rand() * (0.3 - 0.15) + 0.15;
p=0.5;
q=[];q1=[];
c1=[];c2=[];
for pp=1:1
filename=append('AA',int2str(pp),'.mat');
load(filename)
B=A1;
[n m]=size(B);
for i=1:n
for j=1:m
if B(i,j)>0
B(i,j)=1;
else B(i,j)=0;
end
end
end
b=eye(m);
b1=eye(n);
b=eye(m);
b1=eye(n);
for i=1:m
for j=1:m
if i==j
b(i,j)=1;
else
b(i,j)=0.0;
end
end
end
for i=1:n
for j=1:n
if i==j
b1(i,j)=1;
else
b1(i,j)=0.0;
end
end
end
k1=sum(B,1);
k2=sum(B,2);
g=rand() * (1.2 - 0.8) + 0.8;
for ii=1:m
B1(:,ii)=(B(:,ii)./(k1(ii)^p))*g;
end
for ii=1:n
B2(ii,:)=(B(ii,:)./(k2(ii)^p))*g;
end
A= diag(rand(1,n)*(1.1-0.8)+0.8) + (rand(n)-eye(n))*(0.05-0.01) + 0.01*eye(n);
% ts=0:0.05:3;
y0=[];
y0=[rand(m,1); rand(n,1)];
y0=y0';
y0=reshape(y0,[1,m+n]);
p0=y0(1:m);
q0=y0(m+1:m+n);
x=p0; % initial pollinator abundance
y=q0; % initial plant abundance
z=0.0001*rand(m,1)'; %initial norm abundance
k3=0.18;
d=0.5;
m3=0.5;
dA=0.6;
muA=0.0001;
muP=0.0001;
l=0.14;
t=0;
t_max =500; dt=0.01;
m1=t_max/dt;
k=0;
k_max=1;
k_n=10;
dk=(k_max-k)/k_n;
%
for jj=1:k_n+1
t=0.0;
while t<t_max
for j=1:m
% for i=1:n
c1(j)=B1(:,j)'*y';
c1(j)=c1(j)/(1+h*c1(j)); % growth due to mutualism for pollinators
end
for j=1:n
c2(j)=B2(j,:)*x';
c2(j)=c2(j)/(1+h*c2(j)); % growth due to mutualism for plants
end
a1= rand() * (0.35 - 0.05) + 0.05;
a2= rand() * (0.35 - 0.05) + 0.05;
B3=A*(x'.*x); % competition faced by pollinators
B4=A*(y'.*y); % competition faced by plants
x=x+ (a1.*x-diag(B3)'+muA+c1.*x)*dt;
y=y+(a2.*y-k.*y-diag(B4)'+muP+c2.*y)*dt;
t=t+dt;
end
q=[q; k x y ];
k=k+dk;
save(append('data_norm_opt_pol',int2str(pp)),'q','-v7.3');
end
end
figure
% plot(q(:,1),q(:,2:26));
plot(q(:,1),q(:,27:51));
% hold on
% plot(q(:,1),q(:,52:76));
PLease complete the RK4 method to integration;

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by