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 CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by