Help with matlab script

Good day, we were given an assignment to complete. I was able to complete most of it but I was having some trouble approaching the end. I'm not very good using MATLAB and was hoping for some help. The problems I was having were:

1) When I run the program, I get no errors but I'm not seeing the values appear in the workspace.

2) I was trying to perform the least squares method after step 5 using while loops but i'm also having some trouble with that.

3) I wanted to have the program in separate m files like having one main m file, then other supporting m files but I'm also having some trouble with that.

Asking for help because im not fully sure how to complete it.

This is what I have so far:

% compute_drag_script
D = [  0;  0.5;  1.1;  1.7;  2.3;  2.9;  3.5;  4.0;  4.6;  5.2;...
    5.7;  6.3;  6.8;  7.4;  7.9;  8.4;  9.0;  9.5; 10.0; 10.6;...
    11.1; 11.6; 12.1; 12.6; 13.1; 13.6; 14.1; 14.6; 15.1; 15.6;...
    16.1; 16.5; 17.0; 17.5; 18.0; 18.4; 18.9; 19.4]; 
V = [108;  107;  106;  105;  104;  104;  103;  102;  101;  100;...
    99;   99;   98;   97;   97;   96;   95;   94;   94;   93;...
    92;   92;   91;   91;   90;   89;   89;   88;   88;   87;...
    87;   86;   86;   85;   85;   84;   83;   82];
interp_speeds = compute_drag(D,V,0.0095,30);

function interp_speeds = compute_drag(...
    true_distances,true_speeds, k,theta_degrees)
arguments
    true_distances(:,1) {mustBeNonnegative}
    true_speeds(:,1) {mustBeNonnegative}
    k (1,1) {mustBeNonnegative}
    %drag_coeff
    theta_degrees(1,1) {mustBeInRange(theta_degrees,0,180)}
end
%% Input validation
if length(true_distances) ~= length(true_speeds)
    error('Distance and speed vectors must be the same length.');
end
if any(diff(true_distances) < 0)
    error('Distance should never go down!');
end

%% Constants and parameters
g = 9.81; % meters per second-squared
m = 0.45; % kilograms
v0_kph = 108; % kilometers per hour, will convert to meters per second
k = 0.00000001;
ro = 1.293 % kg/(m)^3
d = 0.7/pi % m

%% Conversions
v0 = v0_kph * 1000 / 3600; % kph to meters per second
% We could also convert theta_guess to radians, but there's no reason we
% have to, since Matlab has degree variants of the trig functions.
%% Initial Conditions
x0 = 0;
y0 = 0;
dxdt0 = v0*cosd(theta_degrees);
dydt0 = v0*sind(theta_degrees);
init_vals = [ x0; y0; dxdt0; dydt0];

%% Step 2
% setup for diffeq that will be fed to ode45:
% val(1) x:    x'   (i.e. val(3))
% val(2) y:    y'   (i.e. val(4))
% val(3) x':   (-k|v|/m) * x'
% val(4) y':   -g - (k|v|/m) * y'

all_diffeq = @(t,VALS) [
    VALS(3);
    VALS(4);
    -k * sqrt( VALS(3)^2 + VALS(4)^2 ) * VALS(3) / m;
    -g - ( k*sqrt(VALS(3)^2+VALS(4)^2) * VALS(4) / m);
    ];
% To get an appropriate time range, I'm using the worst case scenario:
% the slowest speed (e.g. 82 kph) taking the ball the full distance
% (e.g. 19.4 meters).
minspeed = min(true_speeds) * 1000 / 3600;
max_necessary_time = (true_distances(end)) / minspeed;
tspan = [0, max_necessary_time];
[t,vals] = ode45(all_diffeq, tspan, init_vals);
x = vals(:,1);
y = vals(:,2);
dxdt = vals(:,3);
dydt = vals(:,4);
% s = sqrt(x.^2 + y.^2);

%% Step 3
s_segments = sqrt(1 + (dydt./dxdt).^2);
s_step = [0;diff(x)];
s = cumsum( s_segments .* s_step );
v = sqrt(dxdt.^2 + dydt.^2);
v_kph = v * 3600 / 1000;

%% Step 4
interp_speeds = interp1(s,v_kph,true_distances);
if any(isnan(interp_speeds))
    warning('NaN encountered. You should increase upper bound on time.');
end
end

% Step 5
while V_interpolated >= 0
    E = cumsum(V_interpolated - V_var).^2
end

% E_theta0 = theta + 0.0001

% E_theta = 

% dE_dtheta = (E_theta - E_base) / 0.0001

% E_k = k + 0.0000008

% dE_dk = (E_k - E_base) / 0.00000001

[Dyuman Joshi]

"When I run the program, I get no errors but I'm not seeing the values appear in the workspace."

Weird, it gives an error here -

% compute_drag_script
D = [  0;  0.5;  1.1;  1.7;  2.3;  2.9;  3.5;  4.0;  4.6;  5.2;...
    5.7;  6.3;  6.8;  7.4;  7.9;  8.4;  9.0;  9.5; 10.0; 10.6;...
    11.1; 11.6; 12.1; 12.6; 13.1; 13.6; 14.1; 14.6; 15.1; 15.6;...
    16.1; 16.5; 17.0; 17.5; 18.0; 18.4; 18.9; 19.4]; 
V = [108;  107;  106;  105;  104;  104;  103;  102;  101;  100;...
    99;   99;   98;   97;   97;   96;   95;   94;   94;   93;...
    92;   92;   91;   91;   90;   89;   89;   88;   88;   87;...
    87;   86;   86;   85;   85;   84;   83;   82];
interp_speeds = compute_drag(D,V,0.0095,30);
function interp_speeds = compute_drag(...
    true_distances,true_speeds, k,theta_degrees)
arguments
  true_distances(:,1) {mustBeNonnegative}
  true_speeds(:,1) {mustBeNonnegative}
  k (1,1) {mustBeNonnegative}
  %drag_coeff
  theta_degrees(1,1) {mustBeInRange(theta_degrees,0,180)}
end
%% Input validation
if length(true_distances) ~= length(true_speeds)
    error('Distance and speed vectors must be the same length.');
end
if any(diff(true_distances) < 0)
    error('Distance should never go down!');
end
%% Constants and parameters
g = 9.81; % meters per second-squared
m = 0.45; % kilograms
v0_kph = 108; % kilometers per hour, will convert to meters per second
k = 0.00000001;
ro = 1.293 % kg/(m)^3
d = 0.7/pi % m
%% Conversions
v0 = v0_kph * 1000 / 3600; % kph to meters per second
% We could also convert theta_guess to radians, but there's no reason we
% have to, since Matlab has degree variants of the trig functions.
%% Initial Conditions
x0 = 0;
y0 = 0;
dxdt0 = v0*cosd(theta_degrees);
dydt0 = v0*sind(theta_degrees);
init_vals = [ x0; y0; dxdt0; dydt0];
%% Step 2
% setup for diffeq that will be fed to ode45:
% val(1) x:    x'   (i.e. val(3))
% val(2) y:    y'   (i.e. val(4))
% val(3) x':   (-k|v|/m) * x'
% val(4) y':   -g - (k|v|/m) * y'
all_diffeq = @(t,VALS) [
    VALS(3);
    VALS(4);
    -k * sqrt( VALS(3)^2 + VALS(4)^2 ) * VALS(3) / m;
    -g - ( k*sqrt(VALS(3)^2+VALS(4)^2) * VALS(4) / m);
    ];
% To get an appropriate time range, I'm using the worst case scenario:
% the slowest speed (e.g. 82 kph) taking the ball the full distance
% (e.g. 19.4 meters).
minspeed = min(true_speeds) * 1000 / 3600;
max_necessary_time = (true_distances(end)) / minspeed;
tspan = [0, max_necessary_time];
[t,vals] = ode45(all_diffeq, tspan, init_vals);
x = vals(:,1);
y = vals(:,2);
dxdt = vals(:,3);
dydt = vals(:,4);
% s = sqrt(x.^2 + y.^2);
%% Step 3
s_segments = sqrt(1 + (dydt./dxdt).^2);
s_step = [0;diff(x)];
s = cumsum( s_segments .* s_step );
v = sqrt(dxdt.^2 + dydt.^2);
v_kph = v * 3600 / 1000;
%% Step 4
interp_speeds = interp1(s,v_kph,true_distances);
if any(isnan(interp_speeds))
    warning('NaN encountered. You should increase upper bound on time.');
end
end
% Step 5
while V_interpolated >= 0

Function definitions in a script must appear at the end of the file.
Move all statements after the "compute_drag" function definition to before the first local function definition.

    E = cumsum(V_interpolated - V_var).^2
end

And as the error states, All function defined in a script must be at the end of the file.

Also, you are using the variables "V_interpolated" and "V_var", but you have not defined them?

"I wanted to have the program in separate m files like having one main m file, then other supporting m files but I'm also having some trouble with that."

Save the function "compute_drag" in the file "compute_drag.m". (Note that a function file must be saved with the name of the first function.) Then open a new file, and insert the values/inputs there and call the function. Save the second file and call it in the workspace.

[Light]

My apologies, I added in step 5 attempting to calculate the values for E. I meant to say when I ran the code up to step 4 there were no errors, my apologies. I was supposed to use interp_speed instead of V_interpolated. I also meant to assign the values for V in the compute_drag function as V_var.

[Light]

Should I just have one other separate file which is the compute_drag.m file or should I have others? Also within the compute_drag file, would it be the function file containing the arguments?

[Light]

I'm also trying to use while loops to perform the least squares method after step 5 but i'm also having a little trouble with that as well.

Answer 1

[chicken vector]

1) Correct me if I am wrong, but with "I'm not seeing the values appear in the workspace", I think you are saying you only see the variables D, V and interp_seeds. To make it simple, you have to imagine functions as black boxes where you feed your input and they return you the corresponding output. Everything happening inside the function is deleted as soon as the output is fed to the running script. If you run the script question1script.m attached to this answer, you will see every variable saved in the workspace. In there I simply run your code as a script instead using a function.

2) To receive more help, you need to explain in more details what do you want to use LSM for. Do you want to find a fit for the propagated data? If so, try to be more specific with your request.

3) Besides a series of paths from which matlab always check for classes, functions, etc., you can select your current folder either manually by clicking in the icon inside the green circle. The selected folder is on the right of the icon:

Or programmatically by typing:

cd C:\Users\Default\AppData

You can save the code of the function compute_drag in a separate .m file and, as long as the file is in the folder in the top bar, Matlab will now from where the code of the function shall be taken.

For more advanced pathing you want to learn which, addpath and genpath, but I suggest you to leave this for the future.

Finally, two remarks:

• You set k as input of the function, but then you overwrite the value inside it, making the input useless. The topions are to either take k out from the list of inputs, or delete the new definition k = 0.00000001; inside the function. If these are intended to be two different things, change the name of one of them.
• Usually you want to define simulation parameters, such as mass and diameter, outside the function and use them as input. This way you function retains more flexibility and it's easier to adjust it to different scenarios.

[Light]

Okay thanks very much. k and theta are supposed to be values I change using the least squares method until I get the desired results which minimizes the error.

I wanted to write a script using the least squares method to find the optimal values of k and theta. I was trying to do that using while loops but I am having some trouble with it.

[Light]

Should I use k as 0.00000001 or k as 0.0095 ?

[chicken vector]

Sorry, I didn't see at first the attached .pdf.

This is a lot to cover for a novice, so I would focus on one thing at a time.

First, a brief comment about step 3.

By using cumsum, you are creating rectangles approximating your curve to compute its integral.

This might be okay, but if you want to increase the precision, have a look at the function trapz, which is most suited for your purpose. In particular, see the example about trapz(X,Y).

Coming to the least square method, let's consider the steps underlined in your .pdf so we can try to understand how they can be translated to code.

1) Start with an initial guess for k and 𝜃

The initial guess is just a metter of choice. Usually, from a shoot you'd exepct an angle of 20-45° and the drag coefficient over one, so we can try to start from:

initialGuess = [30, 1.5];

2) Calculate the error using steps 1-5 above. Call this E_base

We use the values from the initial guess as parameters of the dynamics. One problem with this is that you define the function handle all_diffeq with a fixed value of k and you can't change it.

Try to use a separate file when you define your dynamic function with the syntax:

dydt = dynamics(t,y,parameters)

An example considering a 1DoF mass-spring system would be:

function dydt = dynamicsMassSpring(t,y,k,m)

% y is a 2x1 vector containing position and velocity
% k is the elastic coefficient of the spring
% m is the attached mass

x = y(1);            % Store displacement wrt rest into x
v = y(2);            % Store velocity into v

dxdt = v * t;        % Compute derivative of position
dvdt = -k*x/m;       % Compute derivative of velocity

dydt = [dxdt; dvdt]

end

You can integrate this dynamics by doing:

tspan = [0 10];      % Time interval
y0 = [3 -1];         % Initial conditions
k = -2;              % Spring constant
m = 3;               % Mass

[t,y] = ode45(@(t,y) dynamicsMassSpring(t,y,k,m), tspan, y0);

Then you have to compare the integrated data with "base" data that you called D and V.

To that that you first need to use interp1, and then compute the error.

The Matlab syntax to compute the sum of the square errors is to do:

error = sum( (propagatedData - baseData).^2 );

The notation .^ is very important because it makes the square element-wise along the differnce betwen the two vector values.

now that we know how to get the error for a specific set of initial guess, we want to generalise this by creating another function.

You already know how to do this, you just need to put everything in a function with the right parameters:

function error = estimateError(initialGuess,tspan,D,V,m)

% 1) Extract theta and k from 'initialGuess'
% 2) Integrate this from using these conditions
% 3) Interpolate propagated solution with avilable time data
% 4) Estimate error

end

10) The objective is to find the values of k and 𝜃 that make dE_dtheta and dE_dk BOTH zero. You can use fsolve to accomplish that.

At this point you are ready to use fsolve to find the corresponding initial solution.

The function fsolve accepts a function and an initial guess, so in your main scripts you can define soem input variables such as tspan,D,V,m.

Then you decide an initial guess.

Finally you can find the the solution that minimises the square error by doing:

% Input data: D V m tspan
...
    
% Initial guess [theta, k]:
initialGuess = [30, 1.5];

% Find result:
solution = fsolve(estimateError, initialGuess);

thetaFinal = solution(1);      % Initial angle
kFinal = solution(2);          % Drag coefficient

[Light]

I'm really sorry for the late response, I had a little emergency.

1) Okay so I should use trapz instead of cumsum for more precision? One of the reasons cumsum was used was because not all the values being integrated is scalar. So even though some of the values are not scalar, I could still use the trapz command?

2) You're saying we should use theta as 30 and k be 1.5?

3) Could we go over your step 2 before we move to the step 10 you pointed?

Im sorry for not fully understanding, i'm trying to take it step by step.

[chicken vector]

1) Every value should be scalar, if not, there might be something wrong.

2) The values of the initial guess are arbitrary, the closer you are to the real solution, the faster the convergence.

3) You already did step 2, you have to wrap everything into a function as I showed:

error = estimateError(initialGuess,tspan,D,V,m)

[Light]

Okay, is there a way to convert non scalar values to scalar?

[Light]

I know you''re probably really busy, but I feel like with a little help, we can finish this today

[chicken vector]

You don't want to integrate a scalar value, but a vector. The function trapz works in the same way of cumsum, the difference is that, instead of multiplying by s_step you feed the values of x and y axes as vectors.

These values are s_segments and the vector of corresponding times, where s_segments is a vector of double, which is the "computer" way of say scalar.

Look at the documentaion of trapz for some examples, you have everything that you need in your code.

x = 1 : 50;
y = exp(x/x(end));
area = trapz(x,y); % Compute integral area with trapz!
str = ['Integral area is equal to ' num2str(area)];

figure;
hold on;
plot(x,y,'k','LineWidth',3);
fill([x(1) x x(end)],[0 y 0], 'c')
text(x((end)-x(1))/2, (y(end)-y(1))/2, str)
xlim([x(1)-1,x(end)+1])

Result:

[Light]

So would (s) be equal to trapz(x, y); ?

[chicken vector]

yes but your x is the time and your y is s_segment

[Light]

okay so s = trapz( tspan , s_segments );

[Light]

%% Step 3

s_segments = sqrt(1 + (dydt./dxdt).^2);

s_step = [0;diff(x)];

s = trapz( tspan , s_segments );

v = sqrt(dxdt.^2 + dydt.^2);

v_kph = v * 3600 / 1000;

?

[chicken vector]

Not exactly. tspan contains initial and final time.

What you want is a vector of times that correspond to the values in s_segments.

You have what you need from the output of your ode:

[t,vals] = ode45(all_diffeq, tspan, init_vals);

The variable t containts the time corresponding to the valeus of s_segments, so:

s = trapz(t, s_segments);

[Torsten]

@chicken vector

s = trapz(t, s_segments);

Be careful here. Integration of

s_segments = sqrt(1 + (dydt./dxdt).^2);

must be performed with respect to x, not with respect to t.

[chicken vector]

My mistake! Thank you @Torsten for pointing that out.

s = trapz(vals(:,1), s_segments);

[Light]

Thanks @Torsten and @chicken vector

%% Step 3

s_segments = sqrt(1 + (dydt./dxdt).^2);

s_step = [0;diff(x)];

s = trapz(VALS(:,1), s_segments);

v = sqrt(dxdt.^2 + dydt.^2);

v_kph = v * 3600 / 1000;

[Light]

% compute_drag_script

D = [ 0; 0.5; 1.1; 1.7; 2.3; 2.9; 3.5; 4.0; 4.6; 5.2;...

5.7; 6.3; 6.8; 7.4; 7.9; 8.4; 9.0; 9.5; 10.0; 10.6;...

11.1; 11.6; 12.1; 12.6; 13.1; 13.6; 14.1; 14.6; 15.1; 15.6;...

16.1; 16.5; 17.0; 17.5; 18.0; 18.4; 18.9; 19.4];

V = [108; 107; 106; 105; 104; 104; 103; 102; 101; 100;...

99; 99; 98; 97; 97; 96; 95; 94; 94; 93;...

92; 92; 91; 91; 90; 89; 89; 88; 88; 87;...

87; 86; 86; 85; 85; 84; 83; 82];

interp_speeds = compute_drag(D,V,0.0095,30);

%% Input validation

if length(true_distances) ~= length(true_speeds)

error('Distance and speed vectors must be the same length.');

end

if any(diff(true_distances) < 0)

error('Distance should never go down!');

end

%% Constants and parameters

g = 9.81; % meters per second-squared

m = 0.45; % kilograms

v0_kph = 108; % kilometers per hour, will convert to meters per second

k = 0.00000001;

ro = 1.293 % kg/(m)^3

d = 0.7/pi % m

%% Conversions

v0 = v0_kph * 1000 / 3600; % kph to meters per second

% We could also convert theta_guess to radians, but there's no reason we

% have to, since Matlab has degree variants of the trig functions.

%% Initial Conditions

x0 = 0;

y0 = 0;

dxdt0 = v0*cosd(theta_degrees);

dydt0 = v0*sind(theta_degrees);

init_vals = [ x0; y0; dxdt0; dydt0];

%% Step 2

% setup for diffeq that will be fed to ode45:

% val(1) x: x' (i.e. val(3))

% val(2) y: y' (i.e. val(4))

% val(3) x': (-k|v|/m) * x'

% val(4) y': -g - (k|v|/m) * y'

all_diffeq = @(t,VALS) [

VALS(3);

VALS(4);

-k * sqrt( VALS(3)^2 + VALS(4)^2 ) * VALS(3) / m;

-g - ( k*sqrt(VALS(3)^2+VALS(4)^2) * VALS(4) / m);

];

% To get an appropriate time range, I'm using the worst case scenario:

% the slowest speed (e.g. 82 kph) taking the ball the full distance

% (e.g. 19.4 meters).

minspeed = min(true_speeds) * 1000 / 3600;

max_necessary_time = (true_distances(end)) / minspeed;

tspan = [0, max_necessary_time];

[t,vals] = ode45(all_diffeq, tspan, init_vals);

x = vals(:,1);

y = vals(:,2);

dxdt = vals(:,3);

dydt = vals(:,4);

% s = sqrt(x.^2 + y.^2);

%% Step 3

s_segments = sqrt(1 + (dydt./dxdt).^2);

s_step = [0;diff(x)];

s = trapz(VALS(:,1), s_segments);

v = sqrt(dxdt.^2 + dydt.^2);

v_kph = v * 3600 / 1000;

%% Step 4

interp_speeds = interp1(s,v_kph,true_distances);

if any(isnan(interp_speeds))

warning('NaN encountered. You should increase upper bound on time.');

end

% Step 5

while interp_speeds >= 0

E = cumsum(V_interpolated - V_var).^2

end

% E_theta0 = theta + 0.0001

% E_theta =

% dE_dtheta = (E_theta - E_base) / 0.0001

% E_k = k + 0.0000008

% dE_dk = (E_k - E_base) / 0.00000001

function interp_speeds = compute_drag(...

true_distances,true_speeds, k,theta_degrees)

arguments

true_distances(:,1) {mustBeNonnegative}

true_speeds(:,1) {mustBeNonnegative}

k (1,1) {mustBeNonnegative}

%drag_coeff

theta_degrees(1,1) {mustBeInRange(theta_degrees,0,180)}

end

end

[Light]

should I start with theta as 30 and k as 1.5?

[chicken vector]

Those values are fine as initial guess, but I don't understand why you want to use a while loop.

You don't have to.

Define as the output of the function the quadratic error and your step 5 is basically done.

Read carefully the instructions above and have a solid try before posting one time at a time.

[Light]

okay i'll give it a shot again.

Thats using:

error = estimateError(initialGuess,tspan,D,V,m)

right?

[chicken vector]

You have to create that function but you have already everything in your code.

You function interp_drag does all the necessary steps, except for computing the error.

Re-organise your function such that the input are (initialGuess,tspan,D,V,m), and the output is the error.

Than you can use fsolve to find the solution.

[Light]

Sorry I was still having some trouble with it. Is this any better so far?

% compute_drag_script

D = [ 0; 0.5; 1.1; 1.7; 2.3; 2.9; 3.5; 4.0; 4.6; 5.2;...

5.7; 6.3; 6.8; 7.4; 7.9; 8.4; 9.0; 9.5; 10.0; 10.6;...

11.1; 11.6; 12.1; 12.6; 13.1; 13.6; 14.1; 14.6; 15.1; 15.6;...

16.1; 16.5; 17.0; 17.5; 18.0; 18.4; 18.9; 19.4];

V = [108; 107; 106; 105; 104; 104; 103; 102; 101; 100;...

99; 99; 98; 97; 97; 96; 95; 94; 94; 93;...

92; 92; 91; 91; 90; 89; 89; 88; 88; 87;...

87; 86; 86; 85; 85; 84; 83; 82];

interp_speeds = compute_drag(D,V,0.0095,30);

%% Input validation

if length(true_distances) ~= length(true_speeds)

error('Distance and speed vectors must be the same length.');

end

if any(diff(true_distances) < 0)

error('Distance should never go down!');

end

%% Constants and parameters

g = 9.81; % meters per second-squared

m = 0.45; % kilograms

v0_kph = 108; % kilometers per hour, will convert to meters per second

k = 0.095;

theta = 30

rho = 1.293 % density of air in kg/(m)^3

d = 0.7/pi % m

initialGuess = [30, 0.095];

%% Conversions

v0 = v0_kph * 1000 / 3600; % kph to meters per second

% We could also convert theta_guess to radians, but there's no reason we

% have to, since Matlab has degree variants of the trig functions.

%% Initial Conditions

x0 = 0;

y0 = 0;

dxdt0 = v0*cosd(theta_degrees);

dydt0 = v0*sind(theta_degrees);

init_vals = [ x0; y0; dxdt0; dydt0];

%% Step 2

% setup for diffeq that will be fed to ode45:

% val(1) x: x' (i.e. val(3))

% val(2) y: y' (i.e. val(4))

% val(3) x': (-k|v|/m) * x'

% val(4) y': -g - (k|v|/m) * y'

all_diffeq = @(t,VALS) [

VALS(3);

VALS(4);

-k * sqrt( VALS(3)^2 + VALS(4)^2 ) * VALS(3) / m;

-g - ( k*sqrt(VALS(3)^2+VALS(4)^2) * VALS(4) / m);

];

% To get an appropriate time range, I'm using the worst case scenario:

% the slowest speed (e.g. 82 kph) taking the ball the full distance

% (e.g. 19.4 meters).

minspeed = min(true_speeds) * 1000 / 3600;

max_necessary_time = (true_distances(end)) / minspeed;

tspan = [0, max_necessary_time];

[t,vals] = ode45(all_diffeq, tspan, init_vals);

x = vals(:,1);

y = vals(:,2);

dxdt = vals(:,3);

dydt = vals(:,4);

% s = sqrt(x.^2 + y.^2);

%% Step 3

s_segments = sqrt(1 + (dydt./dxdt).^2);

s_step = [0;diff(x)];

s = trapz(VALS(:,1), s_segments);

v = sqrt(dxdt.^2 + dydt.^2);

v_kph = v * 3600 / 1000;

%% Step 4

interp_speeds = interp1(s,v_kph,true_distances);

if any(isnan(interp_speeds))

warning('NaN encountered. You should increase upper bound on time.');

end

% Step 5

function error = estimateError(initialGuess,tspan,D,V,m)

end

solution = fsolve (estimateError, initialGuess, args= (tspan, D, V, m))

% E_theta0 = theta + 0.0001

% E_theta =

% dE_dtheta = (E_theta - E_base) / 0.0001

% E_k = k + 0.0000008

% dE_dk = (E_k - E_base) / 0.00000001

% Calculation of the drag coefficient

CD = 8*k/(rho*d^2*pi);

function interp_speeds = compute_drag(...

true_distances,true_speeds, k,theta_degrees)

arguments

true_distances(:,1) {mustBeNonnegative}

true_speeds(:,1) {mustBeNonnegative}

k (1,1) {mustBeNonnegative}

%drag_coeff

theta_degrees(1,1) {mustBeInRange(theta_degrees,0,180)}

end

end

[Light]

Im trying really hard to finish the script

[Torsten]

There must remain some correlation between the capabilites of a person doing his/her homework and the points he/she gets for his/her work. So you shouldn't be sad not having mastered the difficult tasks.

[Light]

I know, but im just trying to understand so I could complete the script

Answer 2

[Image Analyst]

You had an extra "end" between step 4 and 5, and V_interpolated is not defined. Was it supposed to be interp_speeds instead of V_interpolated?

% compute_drag_script
D = [  0;  0.5;  1.1;  1.7;  2.3;  2.9;  3.5;  4.0;  4.6;  5.2;...
    5.7;  6.3;  6.8;  7.4;  7.9;  8.4;  9.0;  9.5; 10.0; 10.6;...
    11.1; 11.6; 12.1; 12.6; 13.1; 13.6; 14.1; 14.6; 15.1; 15.6;...
    16.1; 16.5; 17.0; 17.5; 18.0; 18.4; 18.9; 19.4]; 
V = [108;  107;  106;  105;  104;  104;  103;  102;  101;  100;...
    99;   99;   98;   97;   97;   96;   95;   94;   94;   93;...
    92;   92;   91;   91;   90;   89;   89;   88;   88;   87;...
    87;   86;   86;   85;   85;   84;   83;   82];
interp_speeds = compute_drag(D,V,0.0095,30);

%==========================================================================================
% Function definition must occur after the script and end with an "end" statement.
function interp_speeds = compute_drag(...
    true_distances,true_speeds, k,theta_degrees)
arguments
    true_distances(:,1) {mustBeNonnegative}
    true_speeds(:,1) {mustBeNonnegative}
    k (1,1) {mustBeNonnegative}
    %drag_coeff
    theta_degrees(1,1) {mustBeInRange(theta_degrees,0,180)}
end
%% Input validation
if length(true_distances) ~= length(true_speeds)
    error('Distance and speed vectors must be the same length.');
end
if any(diff(true_distances) < 0)
    error('Distance should never go down!');
end
%% Constants and parameters
g = 9.81; % meters per second-squared
m = 0.45; % kilograms
v0_kph = 108; % kilometers per hour, will convert to meters per second
k = 0.00000001;
ro = 1.293 % kg/(m)^3
d = 0.7/pi % m
%% Conversions
v0 = v0_kph * 1000 / 3600; % kph to meters per second
% We could also convert theta_guess to radians, but there's no reason we
% have to, since Matlab has degree variants of the trig functions.
%% Initial Conditions
x0 = 0;
y0 = 0;
dxdt0 = v0*cosd(theta_degrees);
dydt0 = v0*sind(theta_degrees);
init_vals = [ x0; y0; dxdt0; dydt0];
%% Step 2
% setup for diffeq that will be fed to ode45:
% val(1) x:    x'   (i.e. val(3))
% val(2) y:    y'   (i.e. val(4))
% val(3) x':   (-k|v|/m) * x'
% val(4) y':   -g - (k|v|/m) * y'
all_diffeq = @(t,VALS) [
    VALS(3);
    VALS(4);
    -k * sqrt( VALS(3)^2 + VALS(4)^2 ) * VALS(3) / m;
    -g - ( k*sqrt(VALS(3)^2+VALS(4)^2) * VALS(4) / m);
    ];
% To get an appropriate time range, I'm using the worst case scenario:
% the slowest speed (e.g. 82 kph) taking the ball the full distance
% (e.g. 19.4 meters).
minspeed = min(true_speeds) * 1000 / 3600;
max_necessary_time = (true_distances(end)) / minspeed;
tspan = [0, max_necessary_time];
[t,vals] = ode45(all_diffeq, tspan, init_vals);
x = vals(:,1);
y = vals(:,2);
dxdt = vals(:,3);
dydt = vals(:,4);
% s = sqrt(x.^2 + y.^2);
%% Step 3
s_segments = sqrt(1 + (dydt./dxdt).^2);
s_step = [0;diff(x)];
s = cumsum( s_segments .* s_step );
v = sqrt(dxdt.^2 + dydt.^2);
v_kph = v * 3600 / 1000;
%% Step 4
interp_speeds = interp1(s,v_kph,true_distances);
if any(isnan(interp_speeds))
    warning('NaN encountered. You should increase upper bound on time.');
end
% Step 5
while V_interpolated >= 0
    E = cumsum(V_interpolated - V_var).^2
end
% E_theta0 = theta + 0.0001
% E_theta = 
% dE_dtheta = (E_theta - E_base) / 0.0001
% E_k = k + 0.0000008
% dE_dk = (E_k - E_base) / 0.00000001
end

[Light]

Yes it was supposed to be interp_speeds instead of V_interpolated

[Image Analyst]

OK, so does that answer your question/fix your issue?

If this Answer solves your original question, then could you please click the "Accept this answer" link to award the answerer with "reputation points" for their efforts in helping you? They'd appreciate it. Thanks in advance. 🙂 Note: you can only accept one answer (so pick the best one) but you can click the "Vote" icon for as many Answers as you want. Voting for an answer will also award reputation points.

[Light]

Okay, thanks, I would like to accept an answer but the issue I was having isnt completely solved. I was trying to use the least squares method using while loops for step 6-7 on the problem sheet and I was having some trouble with that.

I was also hoping to get a way to have the values show like when using a live script, because I cant fully see the values on the right hand side of the written code right now.

If I could get some help with that, that'd be really great