Manual Runge-Kutta for system of two ODEs.
조회 수: 46 (최근 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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/516862/image.jpeg)
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);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/516882/image.png)
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에서 Numerical Integration and Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!