フィルターのクリア

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

3 ビュー (過去 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 ExchangeNumeric Solvers についてさらに検索

製品


リリース

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by