Manual Runge-Kutta for system of two ODEs.
38 ビュー (過去 30 日間)
古いコメントを表示
I am struggling to obtain the correct graph for the system of ODEs as follows:
x'=-y+6x, y'=-y+4x, between t=0,0.7
I can obtain the correct graph using Euler's method, as seen here:
But cannot do the same for a manual Runge Kutta method. And I don't want to use the integrated ode45 functions if I don't have to. What am I doing wrong? My code is below:
clear,clc
h = 0.1
t_beg = 0
t_end = 0.7
x_initial= 0.5
y_initial= -0.5
F_tx=@(x,y)(-x+6*y);
F_ty=@(x,y)(-y+4*x);
t=t_beg:h:t_end;
x=zeros(1,length(t));
x(1)=x_initial;
y=zeros(1,length(t));
y(1)=y_initial;
for i=1:(length(t)-1)
kx1 = F_tx(t(i),x(i));
kx2 = F_tx(t(i)+0.5*h,x(i)+0.5*h*kx1);
kx3 = F_tx((t(i)+0.5*h),(x(i)+0.5*h*kx2));
kx4 = F_tx((t(i)+h),(x(i)+kx3*h));
x(i+1) = x(i) + (1/6)*(kx1+2*kx2+2*kx3+kx4)*h;
ky1 = F_ty(t(i),y(i));
ky2 = F_ty(t(i)+0.5*h,y(i)+0.5*h*ky1);
ky3 = F_ty((t(i)+0.5*h),(y(i)+0.5*h*ky2));
ky4 = F_ty((t(i)+h),(y(i)+ky3*h));
y(i+1) = y(i) + (1/6)*(ky1+2*ky2+2*ky3+ky4)*h;
end
% plot(x,y)
figure(1)
plot(t,y)
hold on
plot(t,x)
0 件のコメント
採用された回答
Alan Stevens
2021 年 2 月 11 日
編集済み: Alan Stevens
2021 年 2 月 11 日
You need to change the order within the loop to
for i=1:(length(t)-1)
kx1 = F_tx(x(i),y(i));
ky1 = F_ty(x(i),y(i));
kx2 = F_tx(x(i)+0.5*h*kx1,y(i)+0.5*h*ky1);
ky2 = F_ty(x(i)+0.5*h*kx1,y(i)+0.5*h*ky1);
kx3 = F_tx((x(i)+0.5*h*kx2),y(i)+0.5*h*ky2);
ky3 = F_ty((x(i)+0.5*h*kx2),(y(i)+0.5*h*ky2));
kx4 = F_tx((x(i)+kx3*h),y(i)+ky3*h);
ky4 = F_ty((x(i)+kx3*h),(y(i)+ky3*h));
x(i+1) = x(i) + (1/6)*(kx1+2*kx2+2*kx3+kx4)*h;
y(i+1) = y(i) + (1/6)*(ky1+2*ky2+2*ky3+ky4)*h;
end
and note that the first argument is x not t.
4 件のコメント
James Tursa
2021 年 2 月 11 日
For three variables x, y, z you still need to respect the order of the k evaluations. Do kx1, ky1, kz1 first. Then do kx2, ky2, kz2. Etc.
Or to get this same effect use the vector approach that Jan has posted.
その他の回答 (2 件)
Jan
2021 年 2 月 11 日
編集済み: Jan
2023 年 7 月 25 日
The diagram looks, yoike you are integrating:
@(t, y) [-y(1)+6*y(2); -y(2)+4*y(1)]
This code produces an equivalent output:
t0 = 0;
tF = 0.7;
x0 = 0.5;
y0 = -0.5;
[t, y] = ode45(@(t,y) [-y(1)+6*y(2); -y(2)+4*y(1)], ...
[t0, tF], [x0, y0]);
figure;
plot(t,y);
But your function to be integrated is something else:
F_tx = @(x,y) (-x + 6 * y);
F_ty = @(x,y) (-y + 4 * x);
Here the function depends on the 1st input, which is t in my code. This is a confusion of "x/y" versus "t/y", whereby your "y" consists of the components x and y.
[EDITED] A working solution:
F = @(t, y) [-y(1) + 6 * y(2); ...
-y(2) + 4 * y(1)];
y = zeros(2, length(t));
y(:, 1) = [x_initial; y_initial];
for i=1:(length(t)-1)
kx1 = F(t(i), y(:, i));
kx2 = F(t(i) + 0.5 * h, y(:, i) + 0.5 * h * kx1);
kx3 = F(t(i) + 0.5 * h, y(:, i) + 0.5 * h * kx2);
kx4 = F(t(i) + h, y(:, i) + kx3 * h);
y(:, i+1) = y(:, i) + (kx1 + 2 * kx2 + 2 * kx3 + kx4) * h / 6;
end
figure()
plot(t, y)
By the way, compare the readability of the code, which contains spaces around the operators.
khalida
2023 年 7 月 24 日
I am struggling to obtain the correct algorithm for the system of ODEs as follows:
y'=z, z'=-((1+5x)/(2x(x+1))z-y^3+5x+11x^2+0.296x^9+0.666x^10+0.5x^11+0.125x^12), between x=0,0.7
y(0)=0, y'(0)=z(0)=0
i need the matlab code for this system of odes
1 件のコメント
Jan
2023 年 7 月 25 日
Please do not append a new question in the section for answers of another question. Open a new thread instead and delete this message here. Post, what you have tried so far.
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!