# RK 6 method giving large errors

4 views (last 30 days)
Sahil Kumar on 8 Aug 2021
Edited: darova on 11 Aug 2021
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');
Sahil Kumar on 8 Aug 2021

darova on 9 Aug 2021
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 CommentsShowHide 1 older comment
darova on 11 Aug 2021
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 on 9 Aug 2021
Did you compare it with solution from ode45 or ode23s in the matlab?
Sahil Kumar on 9 Aug 2021

R2015a

### Community Treasure Hunt

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

Start Hunting!