How to properly create for loop with if statement

Hi, how do you create a for or while loop with an if statement which adds 360 to an existing variable l, if l is negative?

[Stephen23]

The simple and efficient approach:

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 = mod(L,360)

L = 1×19

180 200 220 240 260 280 300 320 340 0 20 40 60 80 100 120 140 160 180

Answer 1

[Star Strider]

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.

.

[Light]

Okay thank you very much, this helps alot.

[Star Strider]

As always, my pleasure!

[Light]

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]

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]

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]

‘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]

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]

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]

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]

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]

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]

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]

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]

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]

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.

[Light]

Right but I meant how would it be written since atan2d is already assigned to L and we're only adding 360 to each negative value.

Would it still be (using random variable V):

V = -180:20:180

V(V < 0) = V(V<0) + 360

[Star Strider]

It would be:

L = atan2d(Y,X);

L(L<0) = L(L<0) + 360;

If I understand correctly what you want to do.

[Light]

Okay, thank you again so very much.

[Light]

Would it be possible to maybe look over my entire code just to ensure it is up to standard?

[Star Strider]

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]

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]

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

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.

.

[Light]

Okay, it runs without error but the plot lines are a bit off

[Star Strider]

I have no idea how to assess that (although I thought the waypoint in the southern Indian Ocean was interesting).

[Light]

I as well 🤔. I have no idea why that is happening though.

[Star Strider]

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]

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]

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.

[Light]

Okay, this is the link for the calculations guide. https://www.aa.quae.nl/en/reken/grootcirkel.html#fig_3

From 1-6

For the calsulations of X, Y and Z right underneath phi, I used cos and sin in the calculations, should I have used cosd and sind?

[Star Strider]

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.

[Light]

Okay thanks very much

[Star Strider]

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

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]

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]

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]

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]

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.

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.’

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]

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]

‘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.

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]

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]

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]

Okay, thank you again so much for all the help and information you provided.

[Star Strider]

As always, my pleasure!

[Light]

Hi Star, hope all is well. Currently have another code and wanted to find out if it would be possible to request your expertise?

[Star Strider]

It would, although I would have to see it before I comment further.

[Light]

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]

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')

figure
plot(VAR_DATA{:,1}, VAR_DATA{:,2}, '.-')
grid on
xlabel('Distance')
ylabel('Velocity')
axis('equal')

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]

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]

No worries. I’m still catching up from yesterday. I’ll take another look at it and make some suggestions.

[Light]

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]

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 = 

Subs = 

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')

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]

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]

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.

It’s not obvious to me how to set this up, given the provided information.

[Light]

The time range is about 0 to 2 seconds.

[Light]

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]

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]

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]

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.)

[Light]

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.

[Torsten]

g = ...;
k = ...;
m = ...;
[T,Y] = ode45(@(t,y)fun(t,y,g,k,m),[0 20],[1 1 1 1])

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]

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]

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)]

[Light]

Im trying to solve mx2 and my2 using ode45, to get x and y

[Torsten]

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))

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]

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]

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]

Okay thanks very much. Does mag_v need to be changed to normv? Or it can remain as mag_v?

[Light]

Also how can you see the values for X and Y? From the plot seen above

[Light]

I apologise for all the questions, im just trying to understand and, complete the problem in the process.

[Torsten]

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:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

[Light]

Okay, thanks very much. My apologies, I actually completed the course but I forgot some of what was taught.

[Light]

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]

I tried using:

f = @ (x) sqrt( 1 + (y(4) ./ y(3)).^2 )
quad (f, 0, Y(:,1))

But I got an error. Or should the initial conditions be used for y(4) and y(3)?

[Torsten]

Additionally solve the ODE

ds/dt = sqrt(1+(y(4)/y(3))^2) 

with inital condition

s(0) = 0.

Thus your y-vector now becomes 5x1 instead of 4x1 with s being the solution component y(5).

[Light]

This has to be done for each value of X(t) and Y(t) right? To get the various values of distance (s)

[Light]

Sorry just a little confused as to how I should proceed in matlab.

Im trying to convert Y(:, 1) to a scalar variable so I could perform integration using the quad function to get the values of s.

[Light]

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]

With respect to using trapz, use cumtrapz instead to get a vector result rather than a scalar result.

[Light]

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]

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.

[Light]

Okay, I used cumtrapz:

x = [Y(:, 1)]'
y = [Y(:, 2)]'

x1 = [0:0.05:0.7528]

f= sqrt(1 + (y(4) ./ y(3)).^2) 

interp1(Y(:,1), Y(:,2), 0.7)


s = cumtrapz(x, f)

but got an error

[Star Strider]

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.

.

[Light]

I tried changing f to a vector but still, no luck

[Torsten]

@Light

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

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]

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?

.

[Light]

Wow thanks alot, I think that solved the problem.

[Torsten]

@Light

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]

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

[Torsten]

If s is the distance travelled until t=x, the correct formula is

s(x) = integral_{t=0}^{t=x} sqrt((dx/dt)^2 + (dy/dt)^2) dt

not

s(x) = integral_{t=0}^{t=x} sqrt(1 + ((dy/dt)/(dx/dt))^2) dt

[Light]

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]

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.

[Torsten]

Sorry, but your formula for s is simply wrong. Not even the physical unit of your s-definition (which should be meters) is correct for the formula.

The distance travelled is the integral over the magnitude of the velocity.