Numerical differentiation and integration

Hello. Programmming is my weakest point and I have a problem regarding numerical analysis. I was given a problem regarding projectile motion with t and s(t) as data list where s(t) is projectile's trajectory in x-y plane and i need to:

1. write a program that calculates the velocity v = ds/dt of the projectile in m/s with 2nd order accuracy for t = 5,15,…,115 s
2. write a program that calculates the velocity using a second-order polynomial interpolation (either Lagrange or Newtonian) at the points t = 74,75,...95 from the values ​​v (t = 75),v (t = 85), and v (t = 95) found in the previous step. From the values ​​of the velocity at the interpolation points, find the time at which the projectile reaches the highest point of the trajectory, at which the velocity becomes minimum.
3. Using the velocity values ​​from the first step of the task, write a program that calculates the integral from t = 0 s to t = 115 s using Simpson's rule. Estimate the error of Simpson's integration in this case and compare it with the deviation of s(t = 115) that you calculated from the table value. Which error dominates: the integration rule or the arithmetic derivation from which the velocity values ​​used in the integration are derived?

I tried to answer the 1st step but i think it's not the wanted result. My code for this is below (any help with the smallest possible code will be useful, thanx in advance) :

clear all,  clc, close all, format short

% Data list
t = [0,5:10:115]; % in sec
s = [0, 4892.3, 14041.3, 22371.3, 29926.1, 36764.4, 42965.4, 48634.8, 53910.3, 58961.8, 63981.6, 69165.5, 74692.1]; % in m
v0 = 1000; % in m/sec
%-------------------------------------------------------------------------------------------------------------------
N = length(t);

v = zeros(1,N);

for i = 2:N
    v(1,i) = (s(i) - s(i-1)) / (t(i) - t(i-1));
end

v(1) = v0;
disp('Velocity in m/s '), v

Velocity in m/s 

v = 1×13

1.0e+03 * 1.0000 0.9785 0.9149 0.8330 0.7555 0.6838 0.6201 0.5669 0.5275 0.5051 0.5020 0.5184 0.5527

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

subplot(2,1,1)
plot(t,s)
title('s - t')
xlabel('t (s)')
ylabel('s (m)')

subplot(2,1,2)
plot(t,v)
title('v - t')
xlabel('t (s)')
ylabel('v (m/s)')

[Torsten]

You compute the velocity with 1st order accuracy. Do you know which formula to use for 2nd order accuracy ? Do you know how to treat the special cases t = 5 and t = 115 ?

[Left Terry]

Hello. Thanks for the reply. I think i made some progress. My new code regarding the 1st step is shown below (I am not sure if i treated the special cases right):

clear all,  clc, close all, format short

% Data list
t = [0,5:10:115]; % in sec
s = [0, 4892.3, 14041.3, 22371.3, 29926.1, 36764.4, 42965.4, 48634.8, 53910.3, 58961.8, 63981.6, 69165.5, 74692.1]; % in m
v0 = 1000; % in m/sec
%-------------------------------------------------------------------------------------------------------------------
N = length(t);
h = 10;
%-------------------------------------------------------------------------------------------------------------------
% 1st order accuracy
v1 = zeros(1,N-1);
for i = 1:N-1
    v1(1,i) = (s(i+1) - s(i)) / (t(i+1) - t(i));
end
v1(end) = v1(end-1);
disp('Velocity in m/s'), v1

Velocity in m/s

v1 = 1×12

978.4600 914.9000 833.0000 755.4800 683.8300 620.1000 566.9400 527.5500 505.1500 501.9800 518.3900 518.3900

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%-----------------------------------------------------------------------------------------------------
% 2nd order accuracy
v2 = zeros(1,N-1);
for i = 1:N-2
    if i == 1
        v2(1,i) = (s(3) - s(1)) / 15;
    else
    v2(1,i) = (s(i+2) - s(i)) / (2*h);
    end    
end

v2(12) = (s(end) - s(end-1))/10;
disp('Velocity in m/s'), v2

Velocity in m/s

v2 = 1×12

936.0867 873.9500 794.2400 719.6550 651.9650 593.5200 547.2450 516.3500 503.5650 510.1850 535.5250 552.6600

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

subplot(2,1,1)
plot(t,s)
title('s - t')
xlabel('t (s)')
ylabel('s (m)')

subplot(2,1,2)
plot(t(2:end),v1)
title('v - t')
xlabel('t (s)')
ylabel('v1 (m/s)')
hold on
plot(t(2:end),v2)

legend('1st order','2nd order')

Answer 1

[Torsten]

t = [0,5:10:115]; 
s = [0, 4892.3, 14041.3, 22371.3, 29926.1, 36764.4, 42965.4, 48634.8, 53910.3, 58961.8, 63981.6, 69165.5, 74692.1]; % in m
n = numel(t);
%1st order
v1(1) = 1000;
for i = 2:n
    v1(1,i) = (s(i) - s(i-1)) / (t(i) - t(i-1));
end
%2nd order
v2(1) = 1000;
v2(2) = s(1)*(t(2)-t(3))/((t(1)-t(2))*(t(1)-t(3)))+...
        s(2)*((t(2)-t(1))+(t(2)-t(3)))/((t(2)-t(1))*(t(2)-t(3)))+...
        s(3)*(t(2)-t(1))/((t(3)-t(1))*(t(3)-t(2)));
for i = 3:n-1
  v2(i) = (s(i+1)-s(i-1))/(t(i+1)-t(i-1));
end
v2(n) = s(n-2)*(t(n)-t(n-1))/((t(n-2)-t(n-1))*(t(n-2)-t(n)))+...
        s(n-1)*(t(n)-t(n-2))/((t(n-1)-t(n-2))*(t(n-1)-t(n)))+...
        s(n)*((t(n)-t(n-2))+(t(n)-t(n-1)))/((t(n)-t(n-2))*(t(n)-t(n-1)));
subplot(2,1,1)
plot(t,s)
title('s - t')
xlabel('t (s)')
ylabel('s (m)')
subplot(2,1,2)
plot(t,v1)
title('v - t')
xlabel('t (s)')
ylabel('v1 (m/s)')
hold on
plot(t,v2)
legend('1st order','2nd order')

[Left Terry]

Thank you so much for your help. Is there any link to an explanation about how to treat the special cases that you mentioned, as well as how to treat the rest of the problem ?

[Torsten]

Case t = t2:

Differentiate the Lagrange polynomial through (t1,s1),(t2,s2) and (t3,s3) and insert t = t2.

Case t = tn:

Differentiate the Lagrange polynomial through (t_n-2,s_n-2),(t_n-1,s_n-1) and (tn,sn) and insert t = tn.

I leave the rest of the problem to you.

[Left Terry]

Hello. I think i solved the problem. I made a change regarding the special case t = 5 sec. Any way here is the code:

clc, clear all, close all

t = [5:10:115]; 
s = [4892.3, 14041.3, 22371.3, 29926.1, 36764.4, 42965.4, 48634.8, 53910.3, 58961.8, 63981.6, 69165.5, 74692.1]; % in m
N = length(t);

% Ex. 1.1: ds/dt calculation from the data set list with 2nd order accuracy
fprintf('Ex. 1 - 1:\n')

Ex. 1 - 1:

fprintf('\n\tThe velocity ds/dt for t in [5,115] is:\n\n')

	The velocity ds/dt for t in [5,115] is:

v(1) = (2*s(2) - 0.5*s(3) - 1.5*s(1)) / 10;
fprintf('\t\tv(t = %3d)\t=\t%.3f m/sec\n',t(1),v(1))

		v(t =   5)	=	955.850 m/sec

for i = 2:N-1
    v(i) = (s(i+1)-s(i-1))/(t(i+1)-t(i-1));
    fprintf('\t\tv(t = %3d)\t=\t%.3f m/sec\n',t(i),v(i))
end

		v(t =  15)	=	873.950 m/sec
		v(t =  25)	=	794.240 m/sec
		v(t =  35)	=	719.655 m/sec
		v(t =  45)	=	651.965 m/sec
		v(t =  55)	=	593.520 m/sec
		v(t =  65)	=	547.245 m/sec
		v(t =  75)	=	516.350 m/sec
		v(t =  85)	=	503.565 m/sec
		v(t =  95)	=	510.185 m/sec
		v(t = 105)	=	535.525 m/sec

v(N) = s(N-2)*(t(N)-t(N-1))/((t(N-2)-t(N-1))*(t(N-2)-t(N)))+...
        s(N-1)*(t(N)-t(N-2))/((t(N-1)-t(N-2))*(t(N-1)-t(N)))+...
        s(N)*((t(N)-t(N-2))+(t(N)-t(N-1)))/((t(N)-t(N-2))*(t(N)-t(N-1)));
fprintf('\t\tv(t = %3d)\t=\t%.3f m/sec\n',t(N),v(N))

		v(t = 115)	=	569.795 m/sec

%--- Comparing graphically 1st and 2nd order accuracy

    % 1st order accuracy
    v1 = zeros(1,N);
    v1(1) = s(1)/t(1);
    for i = 2:N
        v1(1,i) = (s(i) - s(i-1)) / (t(i) - t(i-1));
    end
   
%     v1(N) = v1(N-1);
    
    plot(t,v1)
    hold on
    plot(t,v)
    title('v-t')
    xlabel('t (sec)')
    ylabel('v (m/sec)')
    legend('1st order','2nd order')
    set(legend,'FontSize',20)

    
% Ex. 1.2: Finding values of v in the domain [75,95] from v(75), v(85) and v(95) using 2nd order interpolation
fprintf('\nEx. 1 - 2:\n')

Ex. 1 - 2:

T = zeros(1,3);
V = zeros(1,3);
Ti = [75:95];

for i = 8:10
    T(1,i-7) = t(i);
    V(1,i-7) = v(i);
end
n = length(T);

V = diag(V)*ones(n,1);
m = length(Ti);
Ti = ones(1,m)*diag(Ti);

L = ones(n,m);

for i = 1:n
    for j = 1:n
        if i ~= j
            L(i,:) = L(i,:).*(Ti - T(j)) / (T(i) - T(j));
        end
    end
end
vi = L'*V;
V = zeros(1,m);

fprintf('\n\tVelocity from t = 75 sec to t = 95 sec is:\n\n')

	Velocity from t = 75 sec to t = 95 sec is:

for i = 1:m
    V(1,i) = vi(i);
    fprintf('\t\tv(t = %d)\t=\t%.3f m/sec\n',Ti(i),V(i));
end

		v(t = 75)	=	516.350 m/sec
		v(t = 76)	=	514.198 m/sec
		v(t = 77)	=	512.241 m/sec
		v(t = 78)	=	510.477 m/sec
		v(t = 79)	=	508.907 m/sec
		v(t = 80)	=	507.532 m/sec
		v(t = 81)	=	506.350 m/sec
		v(t = 82)	=	505.363 m/sec
		v(t = 83)	=	504.570 m/sec
		v(t = 84)	=	503.970 m/sec
		v(t = 85)	=	503.565 m/sec
		v(t = 86)	=	503.354 m/sec
		v(t = 87)	=	503.337 m/sec
		v(t = 88)	=	503.513 m/sec
		v(t = 89)	=	503.884 m/sec
		v(t = 90)	=	504.449 m/sec
		v(t = 91)	=	505.208 m/sec
		v(t = 92)	=	506.161 m/sec
		v(t = 93)	=	507.309 m/sec
		v(t = 94)	=	508.650 m/sec
		v(t = 95)	=	510.185 m/sec

%--- Finding the time where the velocity is minimum
[minV, minVidTi] = min(V);
fprintf('\n\tMinimum velocity occurs at:\n\n\t\tt(v_min) =  %d sec\n',Ti(minVidTi))

	Minimum velocity occurs at:

		t(v_min) =  87 sec

% Ex. 1.3: Integration using Simpson's rule in [0,115]
fprintf('\nEx. 1 - 3:\n')

Ex. 1 - 3:

t = [0,5:10:115];
v = [1000,v];
h1 = 115/2;
h2 = 115/3;
h3 = 115/12;

s1 = (h1/3)*(v(1) + 4*v(7) + v(end));
s2 = (3*h2/8)*(v(1) + 3*v(6) + 3*v(9) + v(end));
s3 = (h3/3)*(v(1) + v(13) + ...
    4*( v(2) + v(4) + v(6) + v(8) + v(10) + v(12) ) + ...
    2*( v(3) + v(5) + v(7) + v(9) + v(11) ));

fprintf('\n\tThe integral of v(t) from t = 0 sec to t = 115 sec is:\n');

	The integral of v(t) from t = 0 sec to t = 115 sec is:

fprintf('\n\t\tBy using Simpson''s 1/3 rule:\n'), fprintf('\n\t\t\ts1 = %.1f m',s1)

		By using Simpson's 1/3 rule:

			s1 = 75590.9 m

fprintf(' with a deviation D1 = %1.1f%% from the given value\n',abs(s1-s(12))/s(12)*100)

 with a deviation D1 = 1.2% from the given value

fprintf('\n\t\tBy using Simpson''s 3/8 rule:\n'), fprintf('\n\t\t\ts2 = %.1f m',s2)

		By using Simpson's 3/8 rule:

			s2 = 72949.4 m

fprintf(' with a deviation D2 = %1.1f%% from the given value\n',abs(s2-s(12))/s(12)*100)

 with a deviation D2 = 2.3% from the given value

fprintf('\n\t\tBy using composite Simpson''s rule:\n'), fprintf('\n\t\t\ts3 = %.1f m',s3)

		By using composite Simpson's rule:

			s3 = 76509.1 m

fprintf(' with a deviation D3 = %1.1f%% from the given value\n',abs(s3-s(12))/s(12)*100)

 with a deviation D3 = 2.4% from the given value

[Torsten]

v1(N) = v1(N-1);

Why don't you use the value for v1(N) computed in the loop before ?

Concerning part (3), I think you are told to use only the composite Simpson's rule. And I can't believe that the result of the composite rule (s3) should be worse than (s1) or (s2). But I didn't check your computation in detail.

Note that the distance between the first three time instants t(1), t(2) and t(3) is 5 and 10, thus not equidistant. You will have to compute the 2nd order polynomial through (t(1),v(1)), (t(2),v(2)) and (t(3),v(3)) and explicitly integrate it for Simpson's 1/3 rule.

Same for Simpson's 3/8 rule which involves (t(1),v(1)), (t(2),v(2)), (t(3),v(3)) and (t(4),v(4)).

But maybe I make things too complicated at the interval end points.

[Left Terry]

You were right, I corrected v1(N). I agree also regarding part 3, s3 should be better than the other two and s2 should be better than s1, but it's the other way around. I have to look through it.

[Left Terry]

I was a missing a term. Now the deviation of s3 is better than before but still greater than the other two.

[Torsten]

t = [0,5:10:115]; 
s = [0, 4892.3, 14041.3, 22371.3, 29926.1, 36764.4, 42965.4, 48634.8, 53910.3, 58961.8, 63981.6, 69165.5, 74692.1]; % in m
n = numel(t);
%1st order
v1(1) = 1000;
for i = 2:n
    v1(1,i) = (s(i) - s(i-1)) / (t(i) - t(i-1));
end
%2nd order
v2(1) = 1000;
v2(2) = s(1)*(t(2)-t(3))/((t(1)-t(2))*(t(1)-t(3)))+...
        s(2)*((t(2)-t(1))+(t(2)-t(3)))/((t(2)-t(1))*(t(2)-t(3)))+...
        s(3)*(t(2)-t(1))/((t(3)-t(1))*(t(3)-t(2)));
for i = 3:n-1
  v2(i) = (s(i+1)-s(i-1))/(t(i+1)-t(i-1));
end
v2(n) = s(n-2)*(t(n)-t(n-1))/((t(n-2)-t(n-1))*(t(n-2)-t(n)))+...
        s(n-1)*(t(n)-t(n-2))/((t(n-1)-t(n-2))*(t(n-1)-t(n)))+...
        s(n)*((t(n)-t(n-2))+(t(n)-t(n-1)))/((t(n)-t(n-2))*(t(n)-t(n-1)));
% Composite Simpson rule
syms tt
p = v2(1)*(tt-t(2))*(tt-t(3))/((t(1)-t(2))*(t(1)-t(3)))+...
    v2(2)*(tt-t(1))*(tt-t(3))/((t(2)-t(1))*(t(2)-t(3)))+...
    v2(3)*(tt-t(1))*(tt-t(2))/((t(3)-t(1))*(t(3)-t(2)));
s0 = double(int(p,tt,t(1),t(3)))

s0 = 1.4047e+04

snum = s0 + 10/3*((v2(3)+4*v2(4)+v2(5))+(v2(5)+4*v2(6)+v2(7))+(v2(7)+4*v2(8)+v2(9))+...
                  (v2(9)+4*v2(10)+v2(11))+(v2(11)+4*v2(12)+v2(13)))

snum = 7.4891e+04

[Left Terry]

Hello. Thank you once again. The new deviation now makes sense. Can you recommend any textbooks regarding the topic of numerical analysis, because the textbooks i have are not very helpful ?

[Torsten]

Maybe this one ?

https://drspmaths.wordpress.com/wp-content/uploads/2019/01/advanced-engineering-mathematics_by-erwin-kreyszig-__10th_edition-1.pdf

What is a "good" book is always a very personal decision. Maybe a Google Search can help:

https://www.google.com/search?client=firefox-b-d&q=recommended+books+on+numerical+analysis

[Left Terry]

Thanks. I will check it out.