RK 6 method giving large errors

Question

Sahil Kumar am 8 Aug. 2021

0
Verknüpfen

Direkter Link zu dieser Frage

https://de.mathworks.com/matlabcentral/answers/894517-rk-6-method-giving-large-errors

Bearbeitet: darova am 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');

3 Kommentare
1 älteren Kommentar anzeigen1 älteren Kommentar ausblenden

Sahil Kumar am 8 Aug. 2021

Bearbeitet: Sahil Kumar am 8 Aug. 2021

The coefficients used were taken from this research paper. An Explicit Sixth-Order Runge-Kutta Formula By H. A. Luther </matlabcentral/answers/uploaded_files/706407/IMG_20210808_210628.jpg>

Sahil Kumar am 8 Aug. 2021

Melden Sie sich an, um zu kommentieren.

Melden Sie sich an, um diese Frage zu beantworten.

Answer 1

darova am 9 Aug. 2021

2
Verknüpfen

Direkter Link zu dieser Antwort

https://de.mathworks.com/matlabcentral/answers/894517-rk-6-method-giving-large-errors#answer_763787

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 Kommentare
Keine anzeigenKeine ausblenden

Sahil Kumar am 9 Aug. 2021

Here, I declared the coefficients separately. Did not do it for the initial few stages as they were just single digits. Also, note that the formula for k1, k2 ... are similar, it is just that they are calculated separately for all the four first order diff equations (fx, fy, fu,fv) and are named respectively (k1x, k1y, k1u,k1v), (k2x, k2y, k2u,k2v), and so on. You may refer to the previous comment for the formulae.

The problem I am facing is that the exact and the numerical solution are different, which should not happen. Is the problem in the code or the logic?

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;
%initial values
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.01; % time step
N=ceil((tf-t(1))/h); % number of iterations
%loop for RK6 calculation
for j=1:N
    t(j+1)=t(j)+h;
    
    x_exact(j+1)=vx*t(j+1); % for exact solution x= vx*t
    y_exact(j+1)=vy*t(j+1)-g*0.5*(t(j+1))^2; %for exact solution y=vy*t + 0.5*g*t^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);
    
    C5 = (7-(21)^0.5)/14; % coefficient for independent var
    %coeff for dependent var 
    C51 = 3*(3*(21)^0.5-7);
    C52 = -8*(7-(21)^0.5);
    C53 = 48*(7-(21)^0.5);
    C54 = -3*(21-(21)^0.5);
    
    k5x=fx(t(j)+C5*h,x(j)+(C51*k1x+C52*k2x+C53*k3x+C54*k4x)/392,u(j)+(C51*k1u+C52*k2u+C53*k3u+C54*k4u)/392);
    k5u=fu(t(j)+C5*h,x(j)+(C51*k1x+C52*k2x+C53*k3x+C54*k4x)/392,u(j)+(C51*k1u+C52*k2u+C53*k3u+C54*k4u)/392);
    k5y=fy(t(j)+C5*h,y(j)+(C51*k1y+C52*k2y+C53*k3y+C54*k4y)/392,v(j)+(C51*k1v+C52*k2v+C53*k3v+C54*k4v)/392);
    k5v=fv(t(j)+C5*h,y(j)+(C51*k1y+C52*k2y+C53*k3y+C54*k4y)/392,v(j)+(C51*k1v+C52*k2v+C53*k3v+C54*k4v)/392);
   
    C6 = (7+(21)^0.5)/14; % coefficient for independent var
    %coeff for dependent var 
    C61 = -5*(231+51*(21)^0.5);
    C62 = -40*(7+(21)^0.5);
    C63 = -320*(21)^0.5;
    C64 = 3*(21+121*(21)^0.5);
    C65 = 392*(6+(21)^0.5);
    
    k6x=fx(t(j)+C6*h,x(j)+(C61*k1x+C62*k2x+C63*k3x+C64*k4x+C65*k5x)/1960,u(j)+(C61*k1u+C62*k2u+C63*k3u+C64*k4u+C65*k5u)/1960);
    k6u=fu(t(j)+C6*h,x(j)+(C61*k1x+C62*k2x+C63*k3x+C64*k4x+C65*k5x)/1960,u(j)+(C61*k1u+C62*k2u+C63*k3u+C64*k4u+C65*k5u)/1960);
    k6y=fy(t(j)+C6*h,y(j)+(C61*k1y+C62*k2y+C63*k3y+C64*k4y+C65*k5y)/1960,v(j)+(C61*k1v+C62*k2v+C63*k3v+C64*k4v+C65*k5v)/1960);
    k6v=fv(t(j)+C6*h,y(j)+(C61*k1y+C62*k2y+C63*k3y+C64*k4y+C65*k5y)/1960,v(j)+(C61*k1v+C62*k2v+C63*k3v+C64*k4v+C65*k5v)/1960);
    
    %coeff for dependent var 
    C71 = 15*(22+7*(21)^0.5);
    C72 = 120;
    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);
    
    k7x=fx(t(j)+h,x(j)+(C71*k1x+C72*k2x+C73*k3x+C74*k4x+C75*k5x+C76*k6x)/180,u(j)+(C71*k1u+C72*k2u+C73*k3u+C74*k4u+C75*k5u+C76*k6u)/180);
    k7u=fu(t(j)+h,x(j)+(C71*k1x+C72*k2x+C73*k3x+C74*k4x+C75*k5x+C76*k6x)/180,u(j)+(C71*k1u+C72*k2u+C73*k3u+C74*k4u+C75*k5u+C76*k6u)/180);
    k7y=fy(t(j)+h,y(j)+(C71*k1y+C72*k2y+C73*k3y+C74*k4y+C75*k5y+C76*k6y)/180,v(j)+(C71*k1v+C72*k2v+C73*k3v+C74*k4v+C75*k5v+C76*k6v)/180);
    k7v=fv(t(j)+h,y(j)+(C71*k1y+C72*k2y+C73*k3y+C74*k4y+C75*k5y+C76*k6y)/180,v(j)+(C71*k1v+C72*k2v+C73*k3v+C74*k4v+C75*k5v+C76*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); % max altitude
rmax = max(x); %range
%plotting results
plot(x,y);
hold on;
plot(x_exact,y_exact);
xlabel('X');
ylabel('Y');

darova am 11 Aug. 2021

Bearbeitet: darova am 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')

Melden Sie sich an, um zu kommentieren.

Answer 2

Wan Ji am 9 Aug. 2021

0
Verknüpfen

Direkter Link zu dieser Antwort

https://de.mathworks.com/matlabcentral/answers/894517-rk-6-method-giving-large-errors#answer_763617

Did you compare it with solution from ode45 or ode23s in the matlab?

1 Kommentar
-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden

Sahil Kumar am 9 Aug. 2021

I don't know about those.

Melden Sie sich an, um zu kommentieren.

RK 6 method giving large errors

3 Kommentare
1 älteren Kommentar anzeigen1 älteren Kommentar ausblenden

Antworten (2)

2 Kommentare
Keine anzeigenKeine ausblenden

1 Kommentar
-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden

Siehe auch

Kategorien

Tags

Produkte

Version

Community Treasure Hunt

RK 6 method giving large errors

3 Kommentare 1 älteren Kommentar anzeigen1 älteren Kommentar ausblenden

Antworten (2)

2 Kommentare Keine anzeigenKeine ausblenden

1 Kommentar -1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden

Siehe auch

Kategorien

Tags

Produkte

Version

Community Treasure Hunt

3 Kommentare
1 älteren Kommentar anzeigen1 älteren Kommentar ausblenden

2 Kommentare
Keine anzeigenKeine ausblenden

1 Kommentar
-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden