Numerical differentiation and integration
226 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Left Terry
am 12 Dez. 2024 um 15:26
Bearbeitet: Torsten
am 14 Dez. 2024 um 21:36
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:
- 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
- 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.
- 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
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)')
2 Kommentare
Torsten
am 12 Dez. 2024 um 19:03
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 ?
Akzeptierte Antwort
Torsten
am 12 Dez. 2024 um 20:18
Bearbeitet: Torsten
am 12 Dez. 2024 um 20:25
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')
Weitere Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!