Runge-kutta results do not align with ODE45 solver, what am I doing wrong?
조회 수: 2 (최근 30일)
이전 댓글 표시
I am working on simulation the movement of three bodies with mass M, position P, velocity V, acceleration A, and force applied in it F.
I am given the initial values of position and velocity to work with, and according to Newton laws I can acquire the force and the right hand side of equation
to implement Runge-Kutta and acquire the position of each body, though when I implement runge-kutta and compare the results to those I get from ODE45, I find them disimilar. What am I doing wrong?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/897355/image.png)
I have attached the file of the instruction I was given and have followed.
The only two files that need to be run is ODE.m and testingrunge.m.
댓글 수: 1
Jan
2022년 2월 17일
편집: Jan
2022년 2월 17일
Just some hints for a simplification:
% Replace
P1=[X(1),X(2),X(3)];
% by
P1 = X(1:3);
% Replace
Y_matrix=cell2mat([{X(4),X(5),X(6)};
{A1(1),A1(2),A1(3)};
{X(10),X(11),X(12)};
{A2(1),A2(2),A2(3)};
{X(16),X(17),X(18)};
{A3(1),A3(2),A3(3)}]);
Y=(Y_matrix(:));
% by
Y = [X(4:6), A1, X(10:12), A2, X(16:18), A3].';
% Replace
F=((-1)*(M*M1).*(P-P1)/((norm(P-P1).^3)))-(1*M*M2.*(P-P2)./((norm(P-P2)).^3));
% by:
F= -M*M1 * (P-P1) / norm(P-P1).^3 - M*M2 * (P-P2) / (norm(P-P2).^3);
채택된 답변
Jan
2022년 2월 17일
편집: Jan
2022년 2월 17일
You are using wrong initial values in the Runge Kutta method:
% X(:, 1) = RHS(0, V0); Nope!
X(:, 1) = V0;
The inital values are the initial values, not the right hand side of them.
I simplified the code to make it easier to debug:
function main
p = 0.05;
P1 = [0;0;0]; V1 = [p;p;p];
P2 = [1;0;0]; V2 = [-p;p;p];
P3 = [0;0;1]; V3 = [-p;-p;-p];
V0 = [P1; V1; P2; V2; P3; V3];
h = 0.01;
t0 = 0;
tEnd = 1;
[t, X] = ode45(@RHS, t0:h:tEnd, V0);
figure;
plot(t, X);
[t2, X2] = rungekutta(h, t0, tEnd, V0);
figure;
plot(t2, X2);
end
function Y = RHS(t, X)
M1 = 1; M2 = 2; M3 = 0.5;
P1 = X(1:3);
P2 = X(7:9);
P3 = X(13:15);
A1 = Force(P1, M1, P2, M2, P3, M3) / M1;
A2 = Force(P2, M2, P1, M1, P3, M3) / M2;
A3 = Force(P3, M3, P1, M1, P2, M2) / M3;
Y = [X(4:6); A1; X(10:12); A2; X(16:18); A3];
end
function F = Force(P, M, P1, M1, P2, M2)
F = -M * M1 * (P-P1) / norm(P-P1).^3 ...
-M * M2 * (P-P2) / norm(P-P2).^3;
end
function [t, X] = rungekutta(h, t0, Tend, V0)
t = t0:h:Tend;
X = zeros(numel(V0), length(t)); % Pre-allocation
X(:, 1) = V0; % Fixed bug here!
for i = 1:length(t)-1
k_1 = h * RHS(t(i), X(:, i));
k_2 = h * RHS(t(i) + 0.5*h, X(:, i) + 0.5 * h * k_1);
k_3 = h * RHS(t(i) + h, X(:, i) - k_1 + 2 * k_2);
X(:, i+1) = X(:, i) + (k_1 + 4 * k_2 + k_3) / 6;
end
end
추가 답변 (0개)
참고 항목
카테고리
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!