How to solve 4th order Runge-Kutta for different initial conditions?
5 ビュー (過去 30 日間)
古いコメントを表示
I have a code that solves the 2 populations with 1 initial conditions and just plot that. How do I make 2D plots for different initial conditions?
%dx/dt = 5*x+1*x*y
%dy/dt = -2*y+1*x*y
clear
clc
f=@(t,z) [5*z(1)+1*z(1)*z(2);
-2*z(2)+1*z(1)*z(2)];
% Initial Conditions
t(1)=0;
z(:,1)=[1,1];
% Step size
s=1;
e=10;
N=e/s;
% Update Loop
for i=1:N
% Update time
t(i+1)=t(i)+s;
% Update for y
k1=f(t(i) ,z(:,i));
k2=f(t(i)+s/2,z(:,i)+s/2*k1);
k3=f(t(i)+s/2,z(:,i)+s/2*k2);
k4=f(t(i)+s ,z(:,i)+s *k3);
z(:,i+1)=z(:,i)+h/6*(k1 + 2*k2 + 2*k3 + k4);
end
1 件のコメント
Davide Masiello
2022 年 4 月 22 日
Your code does not run because h is not defined.
Moreover, what exactly would you like to plot in 2D?
採用された回答
Alan Stevens
2022 年 4 月 22 日
Here's a rather crude method (together with some corrections). You should be able to turn this into a much more elegant version:
%dx/dt = 5*x+1*x*y
%dy/dt = -2*y+1*x*y
f=@(t,z) [5*z(1)+1*z(1)*z(2);
-2*z(2)+1*z(1)*z(2)];
% Example initial conditions
x0 = [1, 0.5, 0.1];
y0 = [-1, -0.5, -0.1];
% Loop through different initial conditions
for j = 1:numel(x0)
% Initial Conditions
z(:,1)=[x0(j),y0(j)]; %%% Make y(0) negative
% Step size
s=0.1; %%% Need much smaller step size
e=10;
N=e/s;
t = zeros(1,N);
% Update Loop
for i=1:N
% Update time
t(i+1)=t(i)+s;
% Update for y
k1=f(t(i) ,z(:,i));
k2=f(t(i)+s/2,z(:,i)+s/2*k1);
k3=f(t(i)+s/2,z(:,i)+s/2*k2);
k4=f(t(i)+s ,z(:,i)+s *k3);
z(:,i+1)=z(:,i)+s/6*(k1 + 2*k2 + 2*k3 + k4); %%% s not h
end
x = z(1,:); y = z(2,:);
figure
subplot(2,1,1)
plot(t,x),grid
xlabel('t'),ylabel('x')
subplot(2,1,2)
plot(t,y),grid
xlabel('t'),ylabel('y')
end
3 件のコメント
Alan Stevens
2022 年 4 月 24 日
Like this?
%dx/dt = 5*x+1*x*y
%dy/dt = -2*y+1*x*y
f=@(t,z) [5*z(1)+1*z(1)*z(2);
-2*z(2)+1*z(1)*z(2)];
% Step size
s=0.1;
e=10;
N=e/s;
% Example initial conditions
x0 = [1, 0.5, 0.1];
y0 = [-1, -0.5, -0.1];
t = zeros(1,N);
x = zeros(numel(x0),N); %%%%%%%%%%
y = zeros(numel(y0),N); %%%%%%%%%%
% Loop through different initial conditions
for j = 1:numel(x0)
z(1,1) = x0(j);
z(2,1) = y0(j);
% Update Loop
for i=1:N-1
% Update time
t(i+1)=t(i)+s;
% Update for y
k1=f(t(i) ,z(:,i));
k2=f(t(i)+s/2,z(:,i)+s/2*k1);
k3=f(t(i)+s/2,z(:,i)+s/2*k2);
k4=f(t(i)+s ,z(:,i)+s *k3);
z(:,i+1)=z(:,i)+s/6*(k1 + 2*k2 + 2*k3 + k4); %%% s not h
end
x(j,:) = z(1,:); y(j,:) = z(2,:); %%%%%%%%%%%%%
end
lbl1 = ['(' num2str(x0(1)) ',' num2str(y0(1)) ')'];
lbl2 = ['(' num2str(x0(2)) ',' num2str(y0(2)) ')'];
lbl3 = ['(' num2str(x0(3)) ',' num2str(y0(3)) ')'];
figure
plot(t,x),grid
xlabel('t'),ylabel('x')
legend(lbl1,lbl2,lbl3)
figure
plot(t,y),grid
xlabel('t'),ylabel('y')
legend(lbl1,lbl2,lbl3)
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!