RK 6 method giving large errors

조회 수: 2(최근 30일)
Sahil Kumar
Sahil Kumar 2021년 8월 8일
편집: darova 2021년 8월 11일
The following code solves for a projectile motion without drag using RK6. I am finding huge variations between values derived from RK6 and that from the exact solution. Can anyone explain why this is happening?
clear all;
%Projectile motion without drag solution using RK6
%x'' = 0, y'' = -g Equations to be solved
g=9.807;
vel=100; th_deg=30; %input
x0=0; y0=0; %initial condition
t0=0; tf=100; %time span
vx=vel*cosd(th_deg); %velocity along x
vy=vel*sind(th_deg); %velocity along y
%transforming second order differential equation to first order
%x'=u represented by fx
%u'=0 represented by fu
%y'=v represented by fy
%v'=-g represented by fv
fx=@(t,x,u) u;
fu=@(t,x,u) 0;
fy=@(t,y,v) v;
fv=@(t,y,v) -g;
t(1)=0;
x(1)=0;y(1)=0;
x_exact(1)=0;y_exact(1)=0; %exact solution
u(1)=vx;v(1)=vy;
h=0.001;
N=ceil((tf-t(1))/h);
for j=1:N
t(j+1)=t(j)+h;
x_exact(j+1)=vx*t(j+1);
y_exact(j+1)=vy*t(j+1)-g*0.5*(t(j+1))^2;
k1x=fx(t(j),x(j),u(j));
k1u=fu(t(j),x(j),u(j));
k1y=fy(t(j),y(j),v(j));
k1v=fv(t(j),y(j),v(j));
k2x=fx(t(j)+h,x(j)+k1x,u(j)+k1u);
k2u=fu(t(j)+h,x(j)+k1x,u(j)+k1u);
k2y=fy(t(j)+h,y(j)+k1y,v(j)+k1v);
k2v=fv(t(j)+h,y(j)+k1y,v(j)+k1v);
k3x=fx(t(j)+h/2,x(j)+((3*k1x+k2x)/8),u(j)+((3*k1u+k2u)/8));
k3u=fu(t(j)+h/2,x(j)+((3*k1x+k2x)/8),u(j)+((3*k1u+k2u)/8));
k3y=fy(t(j)+h/2,y(j)+((3*k1y+k2y)/8),v(j)+((3*k1v+k2v)/8));
k3v=fv(t(j)+h/2,y(j)+((3*k1y+k2y)/8),v(j)+((3*k1v+k2v)/8));
k4x=fx(t(j)+2*h/3,x(j)+(8*k1x+2*k2x+8*k3x)/27,u(j)+(8*k1u+2*k2u+8*k3u)/27);
k4u=fu(t(j)+2*h/3,x(j)+(8*k1x+2*k2x+8*k3x)/27,u(j)+(8*k1u+2*k2u+8*k3u)/27);
k4y=fy(t(j)+2*h/3,y(j)+(8*k1y+2*k2y+8*k3y)/27,v(j)+(8*k1v+2*k2v+8*k3v)/27);
k4v=fv(t(j)+2*h/3,y(j)+(8*k1y+2*k2y+8*k3y)/27,v(j)+(8*k1v+2*k2v+8*k3v)/27);
k5x=fx(t(j)+(7-(21)^0.5)*h/14,x(j)+(3*(3*(21)^0.5-7)*k1x-8*(7-(21)^0.5)*k2x+48*(7-(21)^0.5)*k3x-3*(21-(21)^0.5)*k4x)/392,u(j)+3*((3*(21)^0.5-7)*k1u-8*(7-(21)^0.5)*k2u+48*(7-(21)^0.5)*k3u-3*(21-(21)^0.5)*k4u)/392);
k5u=fu(t(j)+(7-(21)^0.5)*h/14,x(j)+(3*(3*(21)^0.5-7)*k1x-8*(7-(21)^0.5)*k2x+48*(7-(21)^0.5)*k3x-3*(21-(21)^0.5)*k4x)/392,u(j)+3*((3*(21)^0.5-7)*k1u-8*(7-(21)^0.5)*k2u+48*(7-(21)^0.5)*k3u-3*(21-(21)^0.5)*k4u)/392);
k5y=fy(t(j)+(7-(21)^0.5)*h/14,y(j)+(3*(3*(21)^0.5-7)*k1y-8*(7-(21)^0.5)*k2y+48*(7-(21)^0.5)*k3y-3*(21-(21)^0.5)*k4y)/392,v(j)+3*((3*(21)^0.5-7)*k1v-8*(7-(21)^0.5)*k2v+48*(7-(21)^0.5)*k3v-3*(21-(21)^0.5)*k4v)/392);
k5v=fv(t(j)+(7-(21)^0.5)*h/14,y(j)+(3*(3*(21)^0.5-7)*k1y-8*(7-(21)^0.5)*k2y+48*(7-(21)^0.5)*k3y-3*(21-(21)^0.5)*k4y)/392,v(j)+3*((3*(21)^0.5-7)*k1v-8*(7-(21)^0.5)*k2v+48*(7-(21)^0.5)*k3v-3*(21-(21)^0.5)*k4v)/392);
k6x=fx(t(j)+(7+(21)^0.5)*h/14,x(j)+(-5*(231+51*(21)^0.5)*k1x-40*(7+(21)^0.5)*k2x-320*(21)^0.5*k3x+3*(21+121*(21)^0.5)*k4x+392*(6+(21)^0.5)*k5x)/1960,u(j)+(-5*(231+51*(21)^0.5)*k1u-40*(7+(21)^0.5)*k2u-320*(21)^0.5*k3u+3*(21+121*(21)^0.5)*k4u+392*(6+(21)^0.5)*k5u)/1960);
k6u=fu(t(j)+(7+(21)^0.5)*h/14,x(j)+(-5*(231+51*(21)^0.5)*k1x-40*(7+(21)^0.5)*k2x-320*(21)^0.5*k3x+3*(21+121*(21)^0.5)*k4x+392*(6+(21)^0.5)*k5x)/1960,u(j)+(-5*(231+51*(21)^0.5)*k1u-40*(7+(21)^0.5)*k2u-320*(21)^0.5*k3u+3*(21+121*(21)^0.5)*k4u+392*(6+(21)^0.5)*k5u)/1960);
k6y=fy(t(j)+(7+(21)^0.5)*h/14,y(j)+(-5*(231+51*(21)^0.5)*k1y-40*(7+(21)^0.5)*k2y-320*(21)^0.5*k3y+3*(21+121*(21)^0.5)*k4y+392*(6+(21)^0.5)*k5y)/1960,v(j)+(-5*(231+51*(21)^0.5)*k1v-40*(7+(21)^0.5)*k2v-320*(21)^0.5*k3v+3*(21+121*(21)^0.5)*k4v+392*(6+(21)^0.5)*k5v)/1960);
k6v=fv(t(j)+(7+(21)^0.5)*h/14,y(j)+(-5*(231+51*(21)^0.5)*k1y-40*(7+(21)^0.5)*k2y-320*(21)^0.5*k3y+3*(21+121*(21)^0.5)*k4y+392*(6+(21)^0.5)*k5y)/1960,v(j)+(-5*(231+51*(21)^0.5)*k1v-40*(7+(21)^0.5)*k2v-320*(21)^0.5*k3v+3*(21+121*(21)^0.5)*k4v+392*(6+(21)^0.5)*k5v)/1960);
k7x=fx(t(j)+h,x(j)+(15*(22+7*(21)^0.5)*k1x+120*k2x+40*(7*(21)^0.5-5)*k3x-63*(3*(21)^0.5-2)*k4x-14*(49+9*(21)^0.5)*k5x+70*(7-(21)^0.5)*k6x)/180,u(j)+(15*(22+7*(21)^0.5)*k1u+120*k2u+40*(7*(21)^0.5-5)*k3u-63*(3*(21)^0.5-2)*k4u-14*(49+9*(21)^0.5)*k5u+70*(7-(21)^0.5)*k6u)/180);
k7u=fu(t(j)+h,x(j)+(15*(22+7*(21)^0.5)*k1x+120*k2x+40*(7*(21)^0.5-5)*k3x-63*(3*(21)^0.5-2)*k4x-14*(49+9*(21)^0.5)*k5x+70*(7-(21)^0.5)*k6x)/180,u(j)+(15*(22+7*(21)^0.5)*k1u+120*k2u+40*(7*(21)^0.5-5)*k3u-63*(3*(21)^0.5-2)*k4u-14*(49+9*(21)^0.5)*k5u+70*(7-(21)^0.5)*k6u)/180);
k7y=fy(t(j)+h,y(j)+(15*(22+7*(21)^0.5)*k1y+120*k2y+40*(7*(21)^0.5-5)*k3y-63*(3*(21)^0.5-2)*k4y-14*(49+9*(21)^0.5)*k5y+70*(7-(21)^0.5)*k6y)/180,v(j)+(15*(22+7*(21)^0.5)*k1v+120*k2v+40*(7*(21)^0.5-5)*k3v-63*(3*(21)^0.5-2)*k4v-14*(49+9*(21)^0.5)*k5v+70*(7-(21)^0.5)*k6v)/180);
k7v=fv(t(j)+h,y(j)+(15*(22+7*(21)^0.5)*k1y+120*k2y+40*(7*(21)^0.5-5)*k3y-63*(3*(21)^0.5-2)*k4y-14*(49+9*(21)^0.5)*k5y+70*(7-(21)^0.5)*k6y)/180,v(j)+(15*(22+7*(21)^0.5)*k1v+120*k2v+40*(7*(21)^0.5-5)*k3v-63*(3*(21)^0.5-2)*k4v-14*(49+9*(21)^0.5)*k5v+70*(7-(21)^0.5)*k6v)/180);
x(j+1)=x(j)+h*(9*k1x + 64*k3x + 49*k5x + 49*k6x + 9*k7x)/180;
u(j+1)=u(j)+h*(9*k1u + 64*k3u + 49*k5u + 49*k6u + 9*k7u)/180;
y(j+1)=y(j)+h*(9*k1y + 64*k3y + 49*k5y + 49*k6y + 9*k7y)/180;
v(j+1)=v(j)+h*(9*k1v + 64*k3v + 49*k5v + 49*k6v + 9*k7v)/180;
if(y(j+1)<0)
break;
end
end
hmax = max(y);
rmax = max(x);
plot(x,y);
hold on;
plot(x_exact,y_exact);
xlabel('X');
ylabel('Y');
  댓글 수: 3
Sahil Kumar
Sahil Kumar 2021년 8월 8일

댓글을 달려면 로그인하십시오.

답변(2개)

darova
darova 2021년 8월 9일
First of all you should rewrite your code to make more readable. The code is very difficult to to check. Make it as clear as possible
t(j+1)=t(j)+h;
x_exact(j+1)=vx*t(j+1);
y_exact(j+1)=vy*t(j+1)-g*0.5*(t(j+1))^2;
tj = t(j);
xj = x(j);
yj = y(j);
uj = u(j);
vj = v(j);
k1x=fx(tj,xj,uj);
% ... analogically
k2x = fx(tj+h,xj+k1x,uj+k1u);
% ...
C3x = (3*k1x+k2x)/8;
C3u = (3*k1u+k2u)/8;
% ... another C3...
k3x=fx(tj+h/2,xj+C3x,uj+C3u);
% ...
% ...
C71 = 15*(22+7*(21)^0.5);
C73 = 40*(7*(21)^0.5-5);
C74 = 63*(3*(21)^0.5-2);
C75 = 14*(49+9*(21)^0.5);
C76 = 70*(7-(21)^0.5);
dx7 = (C71*k1x+120*k2x+C73*k3x-C74*k4x-C75*k5x+C76*k6x)/180;
du7 = (C71*k1u+120*k2u+C73*k3u-C74*k4u-C75*k5u+C76*k6u)/180;
k7x = fx(tj+h, xj+dx7, uj+du7);
% ...
  댓글 수: 2
darova
darova 2021년 8월 11일
I couldn't fidn a mistake. Here is comparison of ode45 and exact solution. Looks like there is a mistake in the algorithm RK6
g=9.807;
vel=100; th_deg=30; %input
x0=0; y0=0; %initial condition
t0=0; tf=100; %time span
vx=vel*cosd(th_deg); %velocity along x
vy=vel*sind(th_deg); %velocity along y
t = linspace(0,9);
x_exact = vx*t; % for exact solution x= vx*t
y_exact = vy*t-g*0.5*t.^2; %for exact solution y=vy*t + 0.5*g*t^2
f = @(t,f) [f(3);f(4);0;-g];
[t1,f1] = ode45(f,t,[0 0 vx vy]);
plot(x_exact,y_exact,'.r')
line(f1(:,1),f1(:,2))
legend('exact solution','ode45')

댓글을 달려면 로그인하십시오.


Wan Ji
Wan Ji 2021년 8월 9일
Did you compare it with solution from ode45 or ode23s in the matlab?
  댓글 수: 1
Sahil Kumar
Sahil Kumar 2021년 8월 9일
I don't know about those.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by