How to solve 4th order Runge-Kutta for different initial conditions?
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Kenny Gee
am 22 Apr. 2022
Kommentiert: Kenny Gee
am 24 Apr. 2022
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 Kommentar
Davide Masiello
am 22 Apr. 2022
Your code does not run because h is not defined.
Moreover, what exactly would you like to plot in 2D?
Akzeptierte Antwort
Alan Stevens
am 22 Apr. 2022
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 Kommentare
Alan Stevens
am 24 Apr. 2022
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)
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Loops and Conditional Statements finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!