Help with matlab script
Show older comments
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
4 Comments
"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
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
on 24 Apr 2023
Light
on 24 Apr 2023
Light
on 24 Apr 2023
Accepted Answer
chicken vector
on 24 Apr 2023
Edited: chicken vector
on 24 Apr 2023
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.
25 Comments
Light
on 24 Apr 2023
Light
on 24 Apr 2023
chicken vector
on 24 Apr 2023
Edited: chicken vector
on 24 Apr 2023
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.
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
on 26 Apr 2023
chicken vector
on 26 Apr 2023
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
on 26 Apr 2023
chicken vector
on 26 Apr 2023
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
on 26 Apr 2023
chicken vector
on 26 Apr 2023
yes but your x is the time and your y is s_segment
Light
on 26 Apr 2023
Light
on 26 Apr 2023
chicken vector
on 26 Apr 2023
Edited: chicken vector
on 26 Apr 2023
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);
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
on 26 Apr 2023
Edited: Torsten
on 26 Apr 2023
s = trapz(vals(:,1), s_segments);
Light
on 26 Apr 2023
Light
on 26 Apr 2023
Light
on 26 Apr 2023
chicken vector
on 26 Apr 2023
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
on 26 Apr 2023
chicken vector
on 26 Apr 2023
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
on 26 Apr 2023
Light
on 26 Apr 2023
Torsten
on 26 Apr 2023
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
on 26 Apr 2023
More Answers (1)
Image Analyst
on 24 Apr 2023
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
3 Comments
Light
on 24 Apr 2023
Image Analyst
on 24 Apr 2023
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.
Categories
Find more on Loops and Conditional Statements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
