Problem using 4th order Runge-Kutta method
3 次查看(过去 30 天)
显示 更早的评论
I am solving 3 equations using Runge-Kutta method. But the solution doesn't match with exact solution. Please suggest where am I making mistake.
clear;
close all;
clc;
Fx = @(t,x,vx,y,vy,z,vz) [vx; 3*t; vy; t; vz; 7*t] ;
dt=0.1 ;
t=linspace(0, 10, 100);
h=dt;
x=zeros(length(t),1);
vx=zeros(length(t),1);
y=zeros(length(t),1);
vy=zeros(length(t),1);
z=zeros(length(t),1);
vz=zeros(length(t),1);
x(1) = 0 ;
vx(1) = 0 ;
y(1) = 0 ;
vy(1) = 0 ;
z(1) = 0 ;
vz(1) = 0 ;
for i=1: length(t)-1
K1 = Fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i));
K2 = Fx(t(i)+0.5*h, x(i)+0.5*h*K1(1), vx(i)+0.5*h*K1(2), y(i)+0.5*h*K1(3), vy(i)+0.5*h*K1(4), z(i)+0.5*h*K1(5), vz(i)+0.5*h*K1(6));
K3 = Fx(t(i)+0.5*h, x(i)+0.5*h*K2(1), vx(i)+0.5*h*K2(2), y(i)+0.5*h*K2(3), vy(i)+0.5*h*K2(4), z(i)+0.5*h*K2(5), vz(i)+0.5*h*K2(6));
K4 = Fx(t(i)+h, x(i)+h*K3(1), vx(i)+h*K3(2), y(i)+h*K3(3), vy(i)+h*K3(4), z(i)+h*K3(5), vz(i)+h*K3(6));
x(i+1) = x(i) + (1/6)*(K1(1)+2*K2(1)+2*K3(1)+2*K4(1))*h;
vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+2*K4(2))*h;
y(i+1) = y(i) + (1/6)*(K1(3)+2*K2(3)+2*K3(3)+2*K4(3))*h;
vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+2*K4(4))*h;
z(i+1) = z(i) + (1/6)*(K1(5)+2*K2(5)+2*K3(5)+2*K4(5))*h;
vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+2*K4(6))*h;
end
plot(t,x,t,y,t,z,'LineWidth',2)
hold on
plot(t,t.^3/2, t,(1/6)*t.^3,t, (7/6)*t.^3,'LineWidth',2)
xlabel('$t$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
ylabel('$x, y, z$', 'FontSize',20,'FontWeight','bold', 'Interpreter','latex');
set(gca,'FontSize',15);
legend('x','y','z','x_exact','y_exact','z_exact')
1 个评论
James Tursa
2023-3-31
编辑:James Tursa
2023-3-31
First thing I suggest you do is rewrite all of your code in a vectorized fashion. E.g., write your function handle as
Fx = @(t,y) [function of t and y(1), y(2), etc.] ;
Where y is a 6-element vector. I would suggest you reorder your terms as follows:
y(1) = x
y(2) = y
y(3) = z
y(4) = vx
y(5) = vy
y(6) = vz
Then your function handle becomes:
Fx = @(t,y) [y(4:6); 3*t; t; 7*t]
This will greatly simplify your downstream code. It will also allow you to directly use this function handle in a call to ode45( ) so you can compare your RK4 results with the MATLAB ode function results.
回答(1 个)
James Tursa
2023-3-31
You have a factor of 2*K4 that shouldn't be there. It should just be K4. E.g.,
Fx = @(t,x,vx,y,vy,z,vz) [vx; 3*t; vy; t; vz; 7*t] ;
dt=0.01 ;
t=linspace(0, 10, ceil(10/dt));
h=dt;
x=zeros(length(t),1);
vx=zeros(length(t),1);
y=zeros(length(t),1);
vy=zeros(length(t),1);
z=zeros(length(t),1);
vz=zeros(length(t),1);
x(1) = 0 ;
vx(1) = 0 ;
y(1) = 0 ;
vy(1) = 0 ;
z(1) = 0 ;
vz(1) = 0 ;
for i=1: length(t)-1
K1 = Fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i));
K2 = Fx(t(i)+0.5*h, x(i)+0.5*h*K1(1), vx(i)+0.5*h*K1(2), y(i)+0.5*h*K1(3), vy(i)+0.5*h*K1(4), z(i)+0.5*h*K1(5), vz(i)+0.5*h*K1(6));
K3 = Fx(t(i)+0.5*h, x(i)+0.5*h*K2(1), vx(i)+0.5*h*K2(2), y(i)+0.5*h*K2(3), vy(i)+0.5*h*K2(4), z(i)+0.5*h*K2(5), vz(i)+0.5*h*K2(6));
K4 = Fx(t(i)+h, x(i)+h*K3(1), vx(i)+h*K3(2), y(i)+h*K3(3), vy(i)+h*K3(4), z(i)+h*K3(5), vz(i)+h*K3(6));
x(i+1) = x(i) + (1/6)*(K1(1)+2*K2(1)+2*K3(1)+K4(1))*h;
vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+K4(2))*h;
y(i+1) = y(i) + (1/6)*(K1(3)+2*K2(3)+2*K3(3)+K4(3))*h;
vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+K4(4))*h;
z(i+1) = z(i) + (1/6)*(K1(5)+2*K2(5)+2*K3(5)+K4(5))*h;
vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+K4(6))*h;
end
plot(t,x,t,y,t,z,'LineWidth',2)
hold on
plot(t,t.^3/2, t,(1/6)*t.^3,t, (7/6)*t.^3,'LineWidth',2)
xlabel('$t$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
ylabel('$x, y, z$', 'FontSize',20,'FontWeight','bold', 'Interpreter','latex');
set(gca,'FontSize',15);
legend('x','y','z','x_exact','y_exact','z_exact')
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!