You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to properly create for loop with if statement
6 views (last 30 days)
Show older comments
Accepted Answer
Star Strider
on 2 Apr 2023
It is necessary to use the index on ‘L’ on both sides of the equation.
L = -180:20:180
L = 1×19
-180 -160 -140 -120 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180
for l = 1:length(L)
if L(l) < 0
L(l) = L(l) + 360;
end
end
L
L = 1×19
180 200 220 240 260 280 300 320 340 0 20 40 60 80 100 120 140 160 180
A one-line version (avoiding the loop) would be —
L = -180:20:180
L = 1×19
-180 -160 -140 -120 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180
L(L<0) = L(L<0) + 360
L = 1×19
180 200 220 240 260 280 300 320 340 0 20 40 60 80 100 120 140 160 180
The one-line version is more efficient.
.
91 Comments
Light
on 2 Apr 2023
I have another question pertaining to the code. So l is the longitude and 360 must be added to the negative values of l for it to be seen on a plot. I attached the code im trying to finish. Would the one lined version of the code to add 360 be appropriate for the rest of the code? L was already declared as atan2d (Y,X)
Star Strider
on 2 Apr 2023
I am not certain where it would be applied, however if is appropriate to the angle, then yes (with appropriate changes to the variables if necessary).
If you have the Mapping Toolbox, there are several functions, such as wrapTo360 and others that could be applicable here.
Light
on 2 Apr 2023
Edited: Light
on 2 Apr 2023
Okay, no i'll have to download the toolbox. Does wrapTo360 convert the negative longitudes to positive ones?
The code you provided earlier would be implemented in line 51, but I meant, since L was already declared as atan2d(Y, X) could it also be declared as L = -180: 20: 180 ?
Star Strider
on 2 Apr 2023
‘Does wrapTo360 convert the negative longitudes to positive ones?’
It’s easiest to do the experiment since it’s available here —
L = -180:20:180;
Lw = wrapTo360(L)
Lw = 1×19
180 200 220 240 260 280 300 320 340 0 20 40 60 80 100 120 140 160 180
So, yes.
‘... but I meant, since L was already declared as atan2d(Y, X) could it also be declared as L = -180: 20: 180 ?’
If you want to, yes. I would leave it as calculated, however, rather than re-defining it.
.
Light
on 2 Apr 2023
Edited: Light
on 2 Apr 2023
Okay, no I dont want to redefine it, its just that after the calculations are executed for L (which is the longitude) there are some negative values and they wont show up on the plot. But I wanted to have a part of the code which adds 360 to L (logitude) so it would be converted to positive values, thus showing up on the plot.
Star Strider
on 2 Apr 2023
Edited: Star Strider
on 2 Apr 2023
You would need to establish the context there first.
I am not certain what you are asking here. I would not substitute the hard-coded ‘L’ vector for the calculated value, since the atan2 arguments could change. I would just leave it as it is.
EDIT — (2 Apr 2023 at 19:23)
I was not aware what ‘L’ is being used for. If you are using it to plot, it may be necessary to edit any other vectors being plotted as functions of it (if any are) as well. If it is a function of another variable, it could be changed without problems.
Light
on 2 Apr 2023
Edited: Light
on 2 Apr 2023
Okay, L is the longitude being plot, so I am leaving L as is (atan2d (Y,X)). Some of the values of L returned as negative, but since they are negative they wouldn't show on the plot because the values are supposed to be between 0-360. Just wanted to find out how can the code provided earlier (adding 360 to all the values < 0) be implemented in line 51 of the image of the code I attached earlier, so 360 would be added to all the values of L which returned as negative, so all the values of L would be positive.
Light
on 2 Apr 2023
Its just taking the code you provided earlier with the for loop and one lined version (with no loop) and fitting it into line 51 of the code I attached earlier. Im just having some trouble doing that. So the values would show on the plot.
Star Strider
on 2 Apr 2023
That sounds reasonable.
What problems are you encountering?
What are the minimum and maximum angles currently being calculated? For example, would it make sense to add 180 to all of them?
Light
on 2 Apr 2023
I'm not fully sure how to implement the for loop or one lined version (with no loop) into line 51 of the code I provided eralier. Its an assignment im having a little trouble with and 360 must be added to the negative values of L (or longitude).
Star Strider
on 2 Apr 2023
O.K.
If you are supposed to add 360 to the negative values of ‘L’ then the one-line version should work. I would put it just after ‘L’ is calculated in the atan2 call.
Light
on 2 Apr 2023
Okay, but how can it be implemented if L was already declared as atan2d? Would it still be: L = -180:20:180
L(L<0) = L(L<0) + 360 ?
Star Strider
on 2 Apr 2023
I created this ‘L’ vector to test and demonstrate the code:
L = -180:20:180
That is the only reason it exists.
In your code, use the ‘L’ vector created by the atan2 call.
Star Strider
on 2 Apr 2023
It would be:
L = atan2d(Y,X);
L(L<0) = L(L<0) + 360;
If I understand correctly what you want to do.
Light
on 2 Apr 2023
Would it be possible to maybe look over my entire code just to ensure it is up to standard?
Star Strider
on 2 Apr 2023
It would, however I only have images of it thus far. I would need the full text, ideally copied and pasted here rather than attached as .m file that I would then have to download and open.
This appears to be a homework assignment, so I will only comment on what I consider to be errors, or suggestions on improvements, and note them. I will not correct the code. If you want me to run it, I would need any data files as well.
Light
on 2 Apr 2023
Edited: Walter Roberson
on 12 Apr 2023
Okay, I will paste it:
% l= longitude, b= latitude
l1 = -96.81
b1 = 32.78
% translating polar coordinates of the first city (Dallas) to cartesian coordinates
x1 = cosd (-96.81)*cosd (32.78)
y1 = sind (-96.81)*cosd (32.78)
z1 = sind (32.78)
l2 = 14.42
b2 = 50.07
% translating polar coordinates of the second city (Prague) to cartesian coordinates
x2 = cosd (14.42)*cosd (50.07)
y2 = sind (14.42)*cosd (50.07)
z2 = sind (50.07)
% calculating angular distance between the two cities (Dallas and Prague)
psi = acosd ((x1*x2)+(y1*y2)+(z1*z2))
x3 = ((x2 - x1)*cosd(psi))/ sind(psi)
y3 = ((y2 - y1)*cosd(psi))/ sind(psi)
z3 = ((z2 - z1)*cosd(psi))/ sind(psi)
% calculating the distance between the two cities (Dallas and Prague)
r = 6371 Distance = 2*pi*r*(psi/360)
% Plotting the great circle path between the two cities
phi = linspace(0, psi, 100)
X = x1*cos(phi) + x3*sin(phi);
Y = y1*cos(phi) + y3*sin(phi);
Z = z1*cos(phi) + z3*sin(phi);
% converting cartesian coordinates x, y, z to polar coordinates l, b % (longitude and latitude)
b = asind(Z)
l = atan2d (Y,X)
l(l < 0) = l(l < 0) + 360
load('topo.mat', 'topo', 'topomap1');
% load earth topography
whos topo topomap1;
% set coordinates of points
Longitude = 0:15:360; Latitude = -72:6:72;
% set background parameters
contour(0:359, -89:90, topo, [0,0], 'w');
axis equal box on
set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90], 'XTick',... [0: 180: 360], 'YTick', [-90: 90: 90]);
%set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90]);
grid; hold on
plot(l, b, 'y.', 'linewidth', 1);
hold off
Star Strider
on 3 Apr 2023
Edited: Star Strider
on 3 Apr 2023
Running it —
% l= longitude, b= latitude
l1 = -96.81
l1 = -96.8100
b1 = 32.78
b1 = 32.7800
% translating polar coordinates of the first city (Dallas) to cartesian coordinates
x1 = cosd (-96.81)*cosd (32.78)
x1 = -0.0997
y1 = sind (-96.81)*cosd (32.78)
y1 = -0.8348
z1 = sind (32.78)
z1 = 0.5414
l2 = 14.42
l2 = 14.4200
b2 = 50.07
b2 = 50.0700
% translating polar coordinates of the second city (Prague) to cartesian coordinates
x2 = cosd (14.42)*cosd (50.07)
x2 = 0.6216
y2 = sind (14.42)*cosd (50.07)
y2 = 0.1598
z2 = sind (50.07)
z2 = 0.7668
% calculating angular distance between the two cities (Dallas and Prague)
psi = acosd ((x1*x2)+(y1*y2)+(z1*z2))
psi = 77.3049
x3 = ((x2 - x1)*cosd(psi))/ sind(psi)
x3 = 0.1625
y3 = ((y2 - y1)*cosd(psi))/ sind(psi)
y3 = 0.2241
z3 = ((z2 - z1)*cosd(psi))/ sind(psi)
z3 = 0.0508
% calculating the distance between the two cities (Dallas and Prague)
r = 6371;
Distance = 2*pi*r*(psi/360)
Distance = 8.5959e+03
% Plotting the great circle path between the two cities
phi = linspace(0, psi, 100)
phi = 1×100
0 0.7809 1.5617 2.3426 3.1234 3.9043 4.6851 5.4660 6.2469 7.0277 7.8086 8.5894 9.3703 10.1512 10.9320 11.7129 12.4937 13.2746 14.0554 14.8363 15.6172 16.3980 17.1789 17.9597 18.7406 19.5214 20.3023 21.0832 21.8640 22.6449
X = x1*cosd(phi) + x3*sind(phi);
Y = y1*cosd(phi) + y3*sind(phi);
Z = z1*cosd(phi) + z3*sind(phi);
% converting cartesian coordinates x, y, z to polar coordinates l, b % (longitude and latitude)
b = asind(Z)
b = 1×100
32.7800 32.8237 32.8606 32.8907 32.9138 32.9301 32.9395 32.9420 32.9376 32.9263 32.9081 32.8831 32.8511 32.8123 32.7667 32.7142 32.6550 32.5890 32.5162 32.4367 32.3505 32.2576 32.1582 32.0522 31.9396 31.8206 31.6951 31.5632 31.4250 31.2804
l = atan2d (Y,X)
l = 1×100
-96.8100 -96.6843 -96.5576 -96.4298 -96.3008 -96.1706 -96.0392 -95.9063 -95.7721 -95.6363 -95.4990 -95.3600 -95.2193 -95.0768 -94.9323 -94.7859 -94.6374 -94.4867 -94.3337 -94.1783 -94.0204 -93.8599 -93.6966 -93.5304 -93.3612 -93.1888 -93.0131 -92.8339 -92.6511 -92.4644
l(l < 0) = l(l < 0) + 360
l = 1×100
263.1900 263.3157 263.4424 263.5702 263.6992 263.8294 263.9608 264.0937 264.2279 264.3637 264.5010 264.6400 264.7807 264.9232 265.0677 265.2141 265.3626 265.5133 265.6663 265.8217 265.9796 266.1401 266.3034 266.4696 266.6388 266.8112 266.9869 267.1661 267.3489 267.5356
load('topo.mat', 'topo', 'topomap1');
% load earth topography
whos topo topomap1;
Name Size Bytes Class Attributes
topo 180x360 518400 double
topomap1 64x3 1536 double
% set coordinates of points
Longitude = 0:15:360; Latitude = -72:6:72;
figure
% set background parameters
contour(0:359, -89:90, topo, [0,0], 'w');
axis equal
box on
set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90], 'XTick',...
[0: 180: 360], 'YTick', [-90: 90: 90]);
%set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90]);
grid; hold on
plot(l, b, 'y.', 'linewidth', 1);
hold off
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1344054/image.png)
I had to insert some linefeeds to clarify the code lines, however it runs without error. I defer to you to determine if it produces the correct result.
.
Star Strider
on 3 Apr 2023
I have no idea how to assess that (although I thought the waypoint in the southern Indian Ocean was interesting).
Star Strider
on 3 Apr 2023
Takling a closer look, it would seem that perhaps half of the flightpath (for lack of a better term for it) is the mirror image of the other half, with respect to both latitude and lonogitude.
I’m not certain what it’s supposed to be.
Where is it suuposed to start and stop? Is the circumnavigation intentional?
Light
on 3 Apr 2023
Oh well the path is supposed to start at Dallas, Texas and end at Prague, Czech Republic. Its just to calculate the distance between the two cities and plot a path between them.
Star Strider
on 3 Apr 2023
Just out of curiosity, I’ll look up the great circle calculations in the morning and see if I can make it work. I may be able to find the error.
Star Strider
on 3 Apr 2023
Thanks. I pulled up the link, and I’ll go through it in a few minutes.
If you are doing everything in degrees, all the relevant trig functions should be in degrees. Since ‘psi’ is in degrees, the ‘X, Y, Z’ functions should also be in degrees. I changed those and re-ran the code, however while improved, it still ends up in Africa. I’ll look at that in a few minutes.
Star Strider
on 3 Apr 2023
As always, my pleasure!
Here is my analysis, and it appears to provide the correct flight path (that should provide a nice view of tha aurora borealis about midway through the flight for those on the left side of the airplane) —
% % % % % GREAT CIRCLE NAVIGATION
% % % great_circle.m
% % % 03-Mar-2023
% % % Star Strider
b1 = 32.78; % StartLat
l1 = -96.81; % StartLon
b2 = 50.07; % EndLat
l2 = 14.42; % EndLon
r = 6371; % Earth Radius (km)
[x1,y1,z1] = sp2ct(b1,l1);
[x2,y2,z2] = sp2ct(b2,l2);
psi = acosd(x1.*x2+y1.*y2+z1.*z2) % Cities Angular Distance (°)
psi = 77.3049
Distance = 2*pi*r*(psi/360);
fprintf('Distance = %.2f km\n',Distance)
Distance = 8595.92 km
phiv = linspace(0,psi,250).'; % Angular Flight Path
[xv,yv,zv] = angdist(x1,y1,z1,x2,y2,z2,phiv,psi);
[lv,bv] = ct2sp(xv,yv,zv);
lv(lv<0) = lv(lv<0) + 360;
figure
load('topo.mat', 'topo', 'topomap1');
% load earth topography
whos topo topomap1;
Name Size Bytes Class Attributes
topo 180x360 518400 double
topomap1 64x3 1536 double
% set coordinates of points
Longitude = 0:15:360; Latitude = -72:6:72;
% set background parameters
contour(0:359, -89:90, topo, [0,0], 'w');
axis equal
box on
set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90], 'XTick',...
[0: 180: 360], 'YTick', [-90: 90: 90]);
%set(gca, 'Color', [0 0 1], 'XLim', [0 360], 'YLim', [-90 90]);
grid; hold on
plot(lv, bv, 'y.', 'linewidth', 1);
hold off
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1344144/image.png)
function [lon,lat] = ct2sp(x,y,z)
lon = atan2d(y,x);
lat = asind(z);
end
function [x,y,z] = angdist(x1,y1,z1,x2,y2,z2,phi,psi)
% x90 = (x2-x1.*cosd(psi))./sind(psi);
% y90 = (y2-y1.*cosd(psi))./sind(psi);
% z90 = (z2-z1.*cosd(psi))./sind(psi);
% x = x1.*cosd(phi).*x90.*sind(phi);
% y = y1.*cosd(phi).*y90.*sind(phi);
% z = z1.*cosd(phi).*z90.*sind(phi);
x = (x1.*sind(psi-phi)+x2.*sind(phi))./sind(psi);
y = (y1.*sind(psi-phi)+y2.*sind(phi))./sind(psi);
z = (z1.*sind(psi-phi)+z2.*sind(phi))./sind(psi);
end
function [x,y,z] = sp2ct(lat,lon)
x = cosd(lon).*cosd(lat);
y = sind(lon).*cosd(lat);
z = sind(lat);
end
I created functions for the repetitive calculations, and also to simplify the code design. (This is simply a personal preference, and I do not intend it to be criticism of your code.) There are also MATLAB functions cart2sph and sph2cart that will do the conversions. I coded the versions of them here that correspond to the provided great circle documentation. I will let you sort this to determine where the errors may be in your code.
For my part, I now have a great circle .m file if I want to use it in the future!
.
Light
on 3 Apr 2023
Wow this is nothing short of amazing. Really well done. I will look at it and go through it meticulously to try to understand what I did wrong with my code. Could the function section go before the plotting of lv and bv?
Star Strider
on 3 Apr 2023
Thank you!
That was actually fun, although it took me a bit longer to code it than I had expected it to.
The functions here are required to be at the end of the calling code (the order in which the appear is not important, and they can be called anywhere within it and by each other if necessary, just as can other such functions), although anonymous functions (these are not such) must appear in the code before they are used.
The plotting section has to be the last part of the calling code, after all the variables it uses are calculated, and before the function section. Note that the functions are all vectorised, and no (explicit) loops are required.
.
Light
on 4 Apr 2023
Right, I'm still going over the code, but could you explain the 'function' and 'angdist' function? Also I saw phi used at the end, but I didnt see where it was declared.
Star Strider
on 4 Apr 2023
There are three functions: ‘ct2sp’, ‘angdist’, and ‘sp2ct’. I named the first and the last so as not to overshadow similar core MATLAB functions.
All the calculations are from the documentation you provided. The ‘angdist’ function is from
, although I initially coded it as
. When I read down that far, it was preferable because it was obviously more efficient. The earlier code I kept in, although commented-out so that it would not execute.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1344469/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1344474/image.png)
The φ value is actually a vector:
phiv = linspace(0,psi,250).'; % Angular Flight Path
and follows the observation:
‘If φ
, then you are in the first city. If φ=ψ, then you are in the second city.’
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1344479/image.png)
so that extends from the departure airport (0) to the destination airport (ψ), where ψ was defined as:
psi = acosd(x1.*x2+y1.*y2+z1.*z2) % Cities Angular Distance (°)
All the calculations are in degrees.
.
Light
on 4 Apr 2023
Okay,its mainly following the guide given by the calculations guide.
Could you just explain the the 'angdist' function and the 'function' subplot at the end?
Star Strider
on 4 Apr 2023
‘Okay,its mainly following the guide given by the calculations guide.’
Pretty much completely following it.
‘Could you just explain the the 'angdist' function and the 'function' subplot at the end?’
I thought I already explained ‘angdist’ as being
applied to each ‘(x1,x2)’, ‘(y1,y2)’, and ‘(z1,z2)’ pair in turn. It is the same relation applied to each pair. It returns the Cartesian coordinates ‘xv’, ‘yv’, and ‘zv’ that are vectors because φ is the vector ‘phiv’. Those are then transformed into the longitude coordinates ‘lv’, and corresponding latitude coordinates ‘bv’, in the ‘ct2sp’ call immediately after that, both of them being vectors as well, because their arguments are vectors. The ‘lv’ vector is further processed by adding 360 to its negative elements. Those are then plotted on the map.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1344494/image.png)
There is no function subplot. The code after the figure call I copied from your existing code because it loads the map and sets its parameters. (I cannot improve on that, and besides it is consistent with the desired result with respect to displaying the flight path.) The only difference between my code and yours is:
plot(lv, bv, 'y.', 'linewidth', 1);
The three functions at the end of my script are the functions I use in the code.
I just now noticed that you did not post the MATLAB version/release that you are using in either of your threads. My code should run without error in all recent releases. If you are having problems with it because of version/release compatibility issues, I might be able to help with that. The online Run feature here and my offline computers all run R2023a.
.
Light
on 4 Apr 2023
Okay thanks very much. I meant how do you use the function command at the end, because im not familiar with it. Im using the R2022b version.
Star Strider
on 4 Apr 2023
As always, my pleasure!
For that, see the documentation section on Create Functions in Files. This is relatively new (becoming available in the last few years). You will definitely have the ability to use them in R2022b.
I learned some things from answering your Question, so thank you for posting it!
Light
on 12 Apr 2023
Hi Star, hope all is well. Currently have another code and wanted to find out if it would be possible to request your expertise?
Light
on 12 Apr 2023
Edited: Walter Roberson
on 21 Apr 2023
Okay no problem, thanks very much. I know I didnt declare the differential equations properly, but this is all I have right now:
%Differential Equations
% m (d2x / dt^2) = -k|v| (dx/dt)
% m (d2y / dt^2) = -mg - k|v| (dy/dt)
% |v| = sqrt ( (dx/ dt)^2 + (dy/ dt)^2 )
= v0*cos, = v0*sin
% k = air resistance coefficient
% 𝜃 (theta) = initial angle of elevation of the ball when it was kicked
% x and y are respectively the horizontal and vertical spatial coordinates
% t is the time measured from the moment the ball is kicked
% |v| = mag_v = is the magnitude of the velocity (speed, changes with time)
% theta = ?
m = 0.45 % in kg
g = 9.81 % gravitational acceleration in m/s^2
v0 = 108 % in km/h
function dBalldt = Football_F(t, Ball)
m = 0.45 % in kg
g = 9.81 % gravitational acceleration in m/s^2
v0 = 108 % in km/h
% theta = ?
ICs when t = 0
x = 0, y = 0
mag_v = sqrt((v0*cos*theta).^2 + (v0*sin*theta).^2)
% dxdt = @ (t,x)[-k*(mag_v*v0*cos*theta); (-m*g)*-k*(mag_v*v0*sin*theta)];
% m_dxdt = (m*dxdt);
dx/dt = Vx
dy/dt = Vy
dVx= -k*(mag_v*v0*cos*theta)
dVy = (-m*g)*-k*(mag_v*v0*sin*theta)
m_dVx = m*(-k*(mag_v*v0*cos*theta))
m_dVy = m*(-m*g)*-k*(mag_v*v0*sin*theta)
[t,x] = ode45 (m_dxdt, [0,20], [0,0,0]);
return
Star Strider
on 13 Apr 2023
Today turned out to be an unscheduled ‘errand day’ so I was away for several hours.
If I am correct, the objective here is to develop a set of differnetial equations and then estimate parameters using the supplied data. If that is the situation, it would be relatively straightforward to code that. See Parameter Estimation for a System of Differential Equations for an example.
Extracting ‘VAR DATA’ to a table and saving it to a .mat file is straightforward.
VAR_DATA = array2table([0 108
0.5 107
1.1 106
1.7 105
2.3 104
2.9 104
3.5 103
4.0 102
4.6 101
5.2 100
5.7 99
6.3 99
6.8 98
7.4 97
7.9 97
8.4 96
9.0 95
9.5 94
10.0 94
10.6 93
11.1 92
11.6 92
12.1 91
12.6 91
13.1 90
13.6 89
14.1 89
14.6 88
15.1 88
15.6 87
16.1 87
16.5 86
17.0 86
17.5 85
18.0 85
18.4 84
18.9 83
19.4 82], 'VariableNames',{'Distance (m)','Speed (kph)'})
VAR_DATA = 38×2 table
Distance (m) Speed (kph)
____________ ___________
0 108
0.5 107
1.1 106
1.7 105
2.3 104
2.9 104
3.5 103
4 102
4.6 101
5.2 100
5.7 99
6.3 99
6.8 98
7.4 97
7.9 97
8.4 96
save('VAR_DATA.mat','VAR_DATA')
L = size(VAR_DATA,1);
t = linspace(0, L-1, L);
figure
stem3(t, VAR_DATA{:,1}, VAR_DATA{:,2}, 'filled')
grid on
xlabel('t')
ylabel('Distance')
zlabel('Velocity')
axis('equal')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1354128/image.png)
figure
plot(VAR_DATA{:,1}, VAR_DATA{:,2}, '.-')
grid on
xlabel('Distance')
ylabel('Velocity')
axis('equal')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1354133/image.png)
I will address this tomorrow. In the interim, using this information and my prototype code, experiment with coding it. Ideally, it would be preferable if there was also a time vector. in addition to distance and velocity. It it possible to get that, or infer it from the data? For example, are the VAR data sampled at a constant sampling interval or frequency? (I don’t know, and I’ve never encountered this previously, so I’ve never explored it. I just enjoy watching the matches.)
However, it might be possible to use the Symbolic Math Toolbox to solve the differential equations in a closed form , with that solution being the objective function for a nonlinear parameter estimation problem. Try that as well.
.
Light
on 13 Apr 2023
Okay no problem thanks very much. Sorry for the late reply, I think we may have to infer from the data so the time range might be from 0 to 2 seconds. t_range (0: 0.1: 2). I have been going through the code and tried making some changes to the one I posted as well.
Star Strider
on 13 Apr 2023
No worries. I’m still catching up from yesterday. I’ll take another look at it and make some suggestions.
Light
on 19 Apr 2023
Edited: Walter Roberson
on 21 Apr 2023
Hi good day Star, I have been going over the code and I have been having some trouble with it. I'm getting errors which im not sure how to resolve, any idea how the code could be fixed? Im really trying my best
Differential Equations
%Differential Equations
% m (d2x / dt^2) = -k|v| (dx/dt)
% m (d2y / dt^2) = -mg - k|v| (dy/dt)
% |v| = sqrt ( (dx/ dt)^2 + (dy/ dt)^2 )
% |v| = mag_v = sqrt((v0*cos(theta).^2 + (v0*sin(theta).^2)))
% k = air resistance coefficient
% 𝜃 (theta) = initial angle of elevation of the ball when it was kicked
% x and y are respectively the horizontal and vertical spatial coordinates
% t is the time measured from the moment the ball is kicked
% |v| = mag_v = is the magnitude of the velocity (speed, changes with time)
global g
%function dBalldt = Football_F(t, Ball)
%Football_F (1,[1 2 3 4])
m = 0.45 % in kg
g = 9.81 % gravitational acceleration in m/s^2
v0 = 108 % in km/h
mag_v = sqrt((v0*cos(theta).^2 + (v0*sin(theta).^2)))
t_range = 0: 0.1: 2;
theta = 30*pi/180;
% v = 10
% k = ?
Vx = v0*cos(theta)
Vy = v0*sin(theta)
ICs = [0 0 v0*cos(theta) v0*sin(theta)]
% when t = 0
% x = 0, y = 0
% dxdt = @ (t,x)[-k*(mag_v*v0*cos*theta); (-m*g)*-k*(mag_v*v0*sin*theta)];
% m_dxdt = (m*dxdt);
dxdt = Vx
dydt = Vy
%dVx= -k*(mag_v*v0*cos*theta)
%dVy = (-m*g)*-k*(mag_v*v0*sin*theta)
m_dVx = m*(-k*(mag_v*v0*cos*theta))
m_dVy = m*(-m*g)*-k*(mag_v*v0*sin*theta)
[t,xyVxVy] = ode45 (Football_F, t_range, ICs);
%extra arrays
x = xyVxVy (:, 1);
y = xyVxVy (:, 2);
Vx = xyVxVy (:, 3);
Vy = xyVxVy (:, 4);
plot (x,y)
grid
[x y]
y_x5_theoretical = interp1 (x, y, 5);
Star Strider
on 19 Apr 2023
I created and uploaded ‘VAR_DATA.mat’ to make this a bit easier.
I am still somewhat lost on this, and have not given it much thought.
syms k t x(t) y(t) Y T
m = sym(450); % g
g = sym(9.81); % m/s^2
v0 = sym(108); % km/h
v = sqrt(diff(x)^2+diff(y)^2);
Eq1 = m*diff(x,2) == -k*v*diff(x);
Eq2 = m*diff(y,2) == -m*g - k*v*diff(y);
[VF,Subs] = odeToVectorField(Eq1,Eq2)
VF =
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1361373/image.png)
Subs =
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1361378/image.png)
deqfcn = matlabFunction(VF,'Vars',{T,Y,k})
deqfcn = function_handle with value:
@(T,Y,k)[Y(2);k.*sqrt(Y(2).^2+Y(4).^2).*Y(2).*(-1.0./4.5e+2)-9.81e+2./1.0e+2;Y(4);k.*sqrt(Y(2).^2+Y(4).^2).*Y(4).*(-1.0./4.5e+2)]
load('VAR_DATA.mat')
VAR_DATA
VAR_DATA = 38×2 table
Distance (m) Speed (kph)
____________ ___________
0 108
0.5 107
1.1 106
1.7 105
2.3 104
2.9 104
3.5 103
4 102
4.6 101
5.2 100
5.7 99
6.3 99
6.8 98
7.4 97
7.9 97
8.4 96
VAR_MTX = table2array(VAR_DATA);
t = linspace(0, 3, size(VAR_MTX,1)).';
B0 = [1; pi/2; VAR_MTX(1,2)];
B = lsqcurvefit(@objfcn, B0, t, VAR_MTX)
Local minimum possible.
lsqcurvefit stopped because the final change in the sum of squares relative to
its initial value is less than the value of the function tolerance.
B = 3×1
50.5159
1.4735
103.6831
xy = objfcn(B,t)
xy = 38×2
10.0711 103.1928
10.0711 103.1606
10.0711 103.0644
10.0711 102.9057
10.0711 102.6865
10.0711 102.4099
10.0711 102.0791
10.0711 101.6982
10.0711 101.2711
10.0711 100.8022
figure
plot3(VAR_MTX(:,1), VAR_MTX(:,2), t, '.')
hold on
plot3(xy(:,1), xy(:,2), t, '-', 'LineWidth',1.5)
grid
xlabel('Distance')
ylabel('Velocity')
zlabel('t')
legend('Data','Fit', 'Location','best')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1361383/image.png)
function xy = objfcn(b,t)
k = b(1);
theta = b(2);
v0 = b(3);
deqfcn = @(T,Y,k)[Y(2);k.*sqrt(Y(2).^2+Y(4).^2).*Y(2).*(-1.0./4.5e+2)-9.81e+2./1.0e+2;Y(4);k.*sqrt(Y(2).^2+Y(4).^2).*Y(4).*(-1.0./4.5e+2)];
ic = [v0*sin(theta) 0 v0*cos(theta) 0];
[t,y] = ode45(@(t,y)deqfcn(t,y,k), t, ic);
xy = y(:,[3 1]);
end
This is how I usually solve these sorts of problems, however I genuinely do not have a good feeling about what the variables are here. We do not have a time vector (although the exact time vector would be helpful), and I am not certain how to model this. I used the Symbolic Math Toolbox to derive the code for the differential equaations simply because it avoids the algebraic errors that I usually make (and that can be frustrating to detect and correct). Experiment with this code to estimate the parameters and correctly determine what the data are that are being modeled.
Note that you can select the variables that ‘objfcn’ returns by defining what columns the ‘xy’ variable selects from the four columns of ‘y’ returned by the ode45 call. These values will be used by lsqcurvefit to fit the data. The ‘names’ of the variables are in the ‘Subs’ output of the odeToVectorField call.
The earlier problem was straightforward aand relatively easy to solve. This one is not.
.
Light
on 20 Apr 2023
True, thanks very much though. I'm trying to follow along and have a general idea of how it should be solved but I am also somewhat lost on this also.
Star Strider
on 20 Apr 2023
As always, my pleasure!
My problem is that it should be obvious that either x or y is distance and either
or
is velocity, I cannot figure out how they relate here. I’ve been away from problems like this for longer than I’d like to admit, concentrating instead on chemical kinetics, physiological modeling, and signal processing.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1362283/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1362288/image.png)
It’s not obvious to me how to set this up, given the provided information.
Light
on 20 Apr 2023
Edited: Light
on 20 Apr 2023
Do you know how steps 1-6 in the "Steps'' section of the problem is to be carried out, using trapz or quad for number 3?
This is what I have thus far:
% k = air resistance coefficient
% 𝜃 (theta) = initial angle of elevation of the ball when it was kicked
% x and y are respectively the horizontal and vertical spatial coordinates
% t is the time measured from the moment the ball is kicked
% |v| = mag_v = is the magnitude of the velocity (speed, changes with time)
global g
%function dBalldt = Football_F(t, Ball)
%Football_F (1,[1 2 3 4])
m = 0.45 % in kg
g = 9.81 % gravitational acceleration in m/s^2
v0 = 108 % in km/h
mag_v = sqrt((v0*cos(theta).^2 + (v0*sin(theta).^2)))
t = 0: 0.1: 2;
theta = 30*pi/180; %initial guess for theta
k = 3 %initial guess for air resistance in N
% v = 10
Vx = v0*cos(theta)
Vy = v0*sin(theta)
ICs = [0 0 v0*cos(theta) v0*sin(theta)]
% when t = 0
% x = 0, y = 0
% dxdt = @ (t,x)[-k*(mag_v*v0*cos*theta); (-m*g)*-k*(mag_v*v0*sin*theta)];
% m_dxdt = (m*dxdt);
dxdt = Vx
dydt = Vy
%dVx= -k*(mag_v*v0*cos*theta)
%dVy = (-m*g)*-k*(mag_v*v0*sin*theta)
m_dVx = m*(-k*(mag_v*v0*cos*theta))
m_dVy = m*(-m*g)*-k*(mag_v*v0*sin*theta)
[t,xyVxVy] = ode45 (Football_F, t, ICs);
%extra arrays
x = xyVxVy (:, 1);
y = xyVxVy (:, 2);
Vx = xyVxVy (:, 3);
Vy = xyVxVy (:, 4);
plot (x,y)
grid
[x y]
y_x5_theoretical = interp1 (x, y, 5)
Star Strider
on 20 Apr 2023
With respect to 3), the derivatives and ‘x’ vectors would be returned from the ode45 call (in my code, columns 2, 3, and 4). I would then compute the integral as trapz(x, sqrt(1 + dy./dx)).
If I could figure out how the results of integrating the differential equations figured into fitting the data, this would likely be straightforward. I still have no idea how to match those.
Light
on 20 Apr 2023
Okay no problem, maybe the computed values for S could be left out for now. Im still a little confused, sorry. Would the code posted earlier be equivalent to steps 1-3?
Star Strider
on 20 Apr 2023
If I could figure out what ‘x’ and ‘y’ corresponded to in the data, then yes. My code is designed to fit the data and estimate the parameters, however that requires knowing how to fit the integrated differential equation to the data, and that continues to elude me. (I’ve tried various conbinations of x and y to the position data, and various combinations of
and
to the velocity data, as well as various other combinations, and none of them have worked.)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1362383/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1362388/image.png)
Light
on 20 Apr 2023
X and y corresponds to the distances in meters if im not mistaken and the use of the interp1 function in step 4 is to get the interpolated values of speed like in the VAR data.
Light
on 20 Apr 2023
Is this,
% m(d2x / dt2) = -k|v|(dx / dt)
% m(d2y / dt2) = -mg -k|v| (dy / dt)
% |v| = sqrt((v0*cos(theta).^2) + (v0*sin(theta).^2)))
?
Light
on 20 Apr 2023
Sorry, im trying to convert :
% m(d2x / dt2) = -k|v|(dx / dt)
% m(d2y / dt2) = -mg -k|v| (dy / dt)
% |v| = sqrt((v0*cos(theta).^2) + (v0*sin(theta).^2)))
To matlab code
I have this thus far:
% dx/dt = x1
% dx1 / dt = x2
% mx2 = (-k)|v|(x1)
% dy/dt = y1
% dy1 / dt = y2
% my2 = (-mg)(-k)|v|(y1)
m = 0.45 % in kg
g = 9.81 % gravitational acceleration in m/s^2
v0 = 108 % in km/h
mag_v = sqrt((v0*cos(theta).^2) + (v0*sin(theta).^2)))
t = 0: 0.1: 2;
theta = 30*pi/180; %initial guess for theta
k = 3 %initial guess for air resistance in N
ICs = [0 0 v0*cos(theta) v0*sin(theta)]
% v = 10
mx2 = @(t,x) [x(1); x(2); -k*mag_v*x(1)]
my2 = @(t,y) [y(1); y(2); -m*g -k*mag_v*y(2)]
Torsten
on 20 Apr 2023
Edited: Torsten
on 20 Apr 2023
m = 0.45; % in kg
g = 9.81; % gravitational acceleration in m/s^2
v0 = 108; % in km/h
v0 = v0*1000/3600; % in m/s
t = 0: 0.1: 2;
theta = 30*pi/180; %initial guess for theta
k = 3; %initial guess for air resistance in N
ICs = [0 0 v0*cos(theta) v0*sin(theta)];
[T,Y] = ode45(@(t,y)fun(t,y,g,k,m),t,ICs);
plot(T,Y(:,1:2))
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1362463/image.png)
function dy = fun(t,y,g,k,m)
normv = sqrt(y(3)^2+y(4)^2);
dy = zeros(4,1);
dy(1) = y(3);
dy(2) = y(4);
dy(3) = -k*normv*y(3)/m;
dy(4) = (-m*g-k*normv*y(4))/m;
end
Light
on 20 Apr 2023
Thanks very much, I dont mean to be difficult but, could it be written using the same notation like what I had in my last comment with:
mx2 = @(t,x) [x(1); x(2); -k*mag_v*x(1)]
my2 = @(t,y) [y(1); y(2); -m*g -k*mag_v*y(2)]
Torsten
on 20 Apr 2023
Edited: Torsten
on 20 Apr 2023
could it be written using the same notation like what I had in my last comment with:
mx2 = @(t,x) [x(1); x(2); -k*mag_v*x(1)]
my2 = @(t,y) [y(1); y(2); -m*g -k*mag_v*y(2)]
No. The 4 equations must be solved together setting the solution vector to z = [x, y, dx/dt, dy/dt].
And because mag_v has to be updated according to the values of dx/dt and dy/dt, the dz/dt vector should be set in a function, not as a function handle. Just as I did it.
Light
on 20 Apr 2023
Okay thanks very much. Does mag_v need to be changed to normv? Or it can remain as mag_v?
Torsten
on 21 Apr 2023
Edited: Torsten
on 21 Apr 2023
Does mag_v need to be changed to normv? Or it can remain as mag_v?
You can name a variable as you like it.
Also how can you see the values for X and Y? From the plot seen above
X(t) is saved in Y(:,1), Y(t) is saved in Y(:,2) as you can see in the plot command:
plot(T,Y(:,1:2))
To see the values, just add the lines
Y(:,1)
Y(:,2)
It seems it's time that you learn MATLAB fundamentals:
Light
on 21 Apr 2023
Okay, thanks very much. My apologies, I actually completed the course but I forgot some of what was taught.
Light
on 21 Apr 2023
Im following along with the steps posted in the problem and im currently at step 3, so, if I wish to solve the integration problem:
s = ∫ √1 + ( (dy/dt) / (dx/dt) ) .^ 2 dx
with the limits on the integral sign being 0 to x
could it be written similar to:
f = @ (x) sqrt( 1 + (y(4) ./ y(3)).^2 )
quad (f, 0, Y(:,1))
?
Light
on 21 Apr 2023
It doesn't have to be quad, just trying to convert Y(:, 1) to a scalar so I can perform the integration to get the values for the distances.
Star Strider
on 21 Apr 2023
With respect to using trapz, use cumtrapz instead to get a vector result rather than a scalar result.
Light
on 21 Apr 2023
Thanks very much Star. Would it then be:
x1 = [0: 0.1: x]
cumtrapz (x1, y)
% since x has multiple variables and y was already declared?
Star Strider
on 21 Apr 2023
Yes. The vector lengths have to be the same, so rather than the colon operator, linspace would be best, unless the ‘x1’ vector was already defined as something else. If it is, use that vector instead.
Star Strider
on 21 Apr 2023
The ‘f’ value is a scalar, so cumtrapz takes that as the dimension argument and returns the error.
I am sort of lost at this point, however if ‘f’ is supposed to be a vector and if the ‘y’ values you want to use are the integrated results of the ode45 call, consider defining ‘f’ as:
f = sqrt(1 + (y(:,4) ./ y(:,3)).^2);
It will be a column vector as well, since its arguments are column vectors.
The interp1 result needs to be saved to a variable. (It will return the value of ‘Y(:,2)’ corresponding to the interpolated value of ‘Y(:,1)’ that equals 0.7.)
I would also leave the values of ‘x’ and ‘y’ as column vectors for consistency, rather than transposing them. Any row vectors created (for example from linspace) need to be transposed to column vectors as well, again for consistency.
.
Torsten
on 21 Apr 2023
Edited: Torsten
on 21 Apr 2023
Am I really that unclear in what I explain to you ?
m = 0.45; % in kg
g = 9.81; % gravitational acceleration in m/s^2
v0 = 108; % in km/h
v0 = v0*1000/3600; % in m/s
t = 0: 0.01: 2;
theta = 30*pi/180; %initial guess for theta
k = 3; %initial guess for air resistance in N
ICs = [0 0 v0*cos(theta) v0*sin(theta) 0]; % <- Add a 0 at position 5 for s(0)=0
[T,Y] = ode45(@(t,y)fun(t,y,g,k,m),t,ICs);
% Plot s
plot(T,Y(:,5)) % <- Plot s
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1363208/image.png)
function dy = fun(t,y,g,k,m)
normv = sqrt(y(3)^2+y(4)^2);
dy = zeros(4,1);
dy(1) = y(3);
dy(2) = y(4);
dy(3) = -k*normv*y(3)/m;
dy(4) = (-m*g-k*normv*y(4))/m;
dy(5) = sqrt(y(3)^2+y(4)^2); %<- Add the (correct) differential equation for s
end
Star Strider
on 21 Apr 2023
Withiout running the complete code, so creating a random matrix for ‘Y’ —
Y = randn(10,4)
Y = 10×4
1.4637 0.6255 -0.7222 -0.2302
-0.1856 -0.2215 -1.8306 -0.5652
2.4529 -0.2015 0.2981 -0.1843
-1.4724 1.1756 -0.8016 0.7237
-2.4397 0.7014 -0.6971 0.8994
0.5510 -0.6227 0.1481 -1.1258
0.7766 -0.9690 -1.4672 -0.6282
1.5979 -0.7435 0.0185 -1.2906
-0.5662 0.6022 0.0303 -0.0747
-1.0421 0.8021 -1.8383 -1.4224
x = Y(:,1);
y = Y(:,2);
f = sqrt(1 + (Y(:,4) ./ Y(:,3)).^2)
f = 10×1
1.0496
1.0466
1.1757
1.3473
1.6324
7.6682
1.0878
69.8429
2.6567
1.2644
s = cumtrapz(x, f)
s = 10×1
0
-1.7287
1.2030
-3.7486
-5.1898
8.7181
9.7058
38.8334
-39.6132
-40.5462
So perhaps something like this?
.
Torsten
on 21 Apr 2023
As noted in your other post, the equation for s is
ds/dt = sqrt((dx/dt)^2 + (dy/dt)^2)
not
ds/dt = sqrt(1+ ((dy/dt)/(dx/dt))^2).
It already follows from the unit of s (meter) that the second version cannot be correct.
If it were, s had unit "second".
Light
on 21 Apr 2023
the formula for s that was given was s = the integral from 0 to x of sqrt(1+ ((dy/dt)/(dx/dt))^2) dx. I probably made a mistake while typing. It was sqrt(1+ ((dy/dt)/(dx/dt))^2) dx
Light
on 21 Apr 2023
okay no problem, but the formula we were told to use has sqrt(1 + ((dy/dt)/(dx/dt))^2) dx. Step 3 on the problem sheet attached, has s = (integral) sqrt(1+ ((dy/dt)/(dx/dt))^2) dx
Light
on 21 Apr 2023
mag_v or |v| (magnitude of v) is given as: sqrt((v0*cos(theta).^2) + (v0*sin(theta).^2))
But the formula for the distance using x and y is: s = (integral) sqrt(1+ ((dy/dt)/(dx/dt))^2) dx
sorry for the confusion. The first is the magnitude of v while the second is the distances based on the values of x and y.
More Answers (0)
See Also
Categories
Find more on Performance 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!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)