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

2 次查看（过去 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 中查找有关 Oceanography and Hydrology 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by