How can I solve a polynomial where coefficients are vector arrays?

Hi there!

I have a system of equations which 4 of them feed into a final equation (code shown below). The first four equations all work beautifully, no issues, even while using large arrays as one of my variables.

The variable I am concerned with here is "psidot" - when I make psidot a single value, say 90, the entire equation works beautiful and my roots come out appropriately and I get the correct answer. However, as soon as I make "psidot" an array, say [90 80 70], my roots get all messed up and I don't get the right answers anymore. I get a lot of negative numbers and certainly none of the right answers.

Essentially what I'm looking for is to evaluate the roots of the polynomial at each of these values of "psidot", individually. I'm just not sure what i'm doing wrong. Can anyone provide some insight?

Thanks so much in advance! Caitlin

v1=674; %ft/sec velocity of aircraft 1

v2=760; %ft/sec velocity of aircraft 2

theta1=0; %degrees reckoning of aircraft 1

theta2=180; %degres reckoning of aircraft 2

psi1=[90 80 70]; %degrees degrees of error (should be array)

psidot=3; %degrees/sec rate of error change

syms t; %time variable

x1=6049; %x position of aircraft 1

y1=12089; %y position of aircraft 1

x2=90735; %x position of aircraft 2

y2=120980; %y position of aircraft 2

d=500; %distance they should be apart for a problem, in feet

k1=v1*cos(theta1+psi1)-v2*cos(theta2); 

k2=((v1/psidot)*(sin(theta1+psi1)-sin(theta1)))-cos(theta1+psi1)*v1.*(psi1/psidot)+x1-x2;

k3=v1*sin(theta1+psi1)-v2*sin(theta2);

k4=(v1/psidot)*(cos(theta1)-cos(theta1+psi1))-sin(theta1+psi1)*v1.*(psi1/psidot)+y1-y2;

k5=k1.^2+k3.^2; %"a" in the polynomial

k6=2*(k1.*k2+k3.*k4); %"b" in the polynomial

k7=k3.^2+k4.^2-500; %"c" in the polynomial

p = [k5 k6 k7]; % creating the polynomial

r = roots(p); %solving for roots

r=real(r)/60 %getting real answers

[Andy L]

Star Strider is correct in his answer below. However I noticed that in k2 and k4 you divide two vectors without using the ./ operator. This means MATLAB will be trying to divide the two matrices as opposed to their elements, which is not possible as the inner dimensions will disagree - i'm assuming this is your error?

[Caitlin]

Thanks - I fixed it :-)

Answer 1

[Star Strider]

If psidot is a vector, use a loop:

psidotvct = [90 80 70]; %degrees/sec rate of error change
for k1 = 1:size(psidotvct,2)
    psidot = psidotvct(k1);

      .... CODE ...

    rr(k1,:) = real(r)/60 %getting real answers
end

This creates the rr array to store one row for every value of psidot. (I changed your code to create rr because it would otherwise cause problems with r just before it in subsequent loop iterations.)

I also don’t understand the ‘syms t’ call. It doesn’t look like you’re using it anywhere. You might consider deleting that line to speed things up.

[Caitlin]

Hi!

Thanks for the help!

Do i need to do a loop for each equation? k1, k2, etc?

Here's what it looks like now and i'm getting an error in my "rr" statement. Also, i realized i made a mistake in my description..i want "psi1" to be an array, not "psidot" psidot is fixed...maybe that makes a difference...sorry about that.

v1=674; %ft/sec velocity of aircraft 1

v2=760; %ft/sec velocity of aircraft 2

theta1=0; %degrees reckoning of aircraft 1

theta2=180; %degres reckoning of aircraft 2

psivct=[90 80 70];%degrees degrees of error (should be array)

psidot=3; %degrees/sec rate of error change

x1=6049; %x position of aircraft 1

y1=12089; %y position of aircraft 1

x2=90735; %x position of aircraft 2

y2=120980; %y position of aircraft 2

d=500; %distance they should be apart for a problem, in feet

for k1 = 1:size(psivct,2);

    psi1=psivct(k1);

k1=v1*cos(theta1+psi1)-v2*cos(theta2);

k2=((v1/psidot)*(sin(theta1+psi1)-sin(theta1)))-cos(theta1+psi1)*v1.*(psi1./psidot)+x1-x2;

k3=v1*sin(theta1+psi1)-v2*sin(theta2);

k4=(v1/psidot)*(cos(theta1)-cos(theta1+psi1))-sin(theta1+psi1)*v1.*(psi1./psidot)+y1-y2;

k5=k1.^2+k3.^2; %"a" in the polynomial

k6=2*(k1.*k2+k3.*k4); %"b" in the polynomial

k7=k3.^2+k4.^2-500; %"c" in the polynomial

p = [k5 k6 k7]; % creating the polynomial

r = roots(p); %solving for roots

rr(k1,:)=real(r)/60 %getting real answers

end;

[Caitlin]

THANK YOU. life saver :-)

[Star Strider]

My pleasure!

I don’t entirely understand what you want to save within your loop. My impression was that you wanted to store the roots given by r for each value of psi1. You don’t need a loop for each k1, k2, etc., but if you want to store them in each iteration of the loop, and they are scalars (as I assume), give them subscripts: k1(c1), k2(c2), ... If they’re row vectors in each iteration and have constant, equal lengths over all the iterations, save them as k1(c1,:)=, if they’re column vectors (same conditions), save them as k1(:,c1)=, etc.

Revisions in my posted code to accommodate your clarification and the fact that I didn’t note that you are using my usual loop counter, k1 as a variable, and the fact that I forgot that roots produces a column vector would be:

psi1vct = [90 80 70]; %degrees degrees of error (should be array)
for c1 = 1:size(psi1vct,2)
    psi1 = psi1vct(c1);
    ... CODE ...
    k1(c1) = ... etc. if they are scalars and you want to save them 
    rr(:,c1) = real(r)/60 %getting real answers
end

If you want to save k1 ... k7 and subscript them, you will have to change all their references in all statements that use them to include their subscripts as well. For instance, your p assignment changes to:

p = [k5(c1) k6(c1) k7(c1)]; % creating the polynomial

If you don’t want to save them from iteration-to-iteration as vectors, the rest of your code, including the p assignment, doesn’t change.

(Noticed later you’re doing research on human factors in aviation, an interest of mine and at one time an area of research.)

[Star Strider]

My pleasure!

[Caitlin]

wow - this is awesome.

i do do human factors in aviation, particularly air traffic control. it's pretty cool :-)

one thing i am working on with this now is i eventually want to work through so that this entire thing is going to be essentially vectors. I have MASSIVE amounts of data dumping in from an FAA data base where i'm looking to evaluate pairs of aircrafts in a sector... so every variable (v1, v2, theta1, theta2, x1, y1, x2, y2) will have one value per aircraft, and there are anywhere from 12-36 aircrafts in the sector, and then i'll be looking at each aircraft for 40 minutes (1 value of each per minute), so it's a lot of data.

To make my other values vectors, do i just add on to my rr statement? ex: rr(:,c1, c2, c3) %if c2 and c3 were my vectors for v1 and v2

and then just add for loops up top to create those arrays? Is that correct?

Thanks again for the awesome help - it's really making a difference for me!

[Star Strider]

My pleasure!

Fun for me, too.

1. One thing I overlooked yesterday but picked up on today is:

theta1=0; %degrees reckoning of aircraft 1
theta2=180; %degres reckoning of aircraft 2

then later you use (for instance):

k1=v1*cos(theta1+psi1)-v2*cos(theta2);

Use cosd (and sind, tand, ... etc.) instead. The cos, sin, and others use radian angles. The cosd, sind, and others use — and their inverses return — degree angles. If you’re using degrees, change your code to use the appropriate functions.

2. I need to understand what k1 ... k7 are and what r represents to answer your question about your other values vectors. (Also, what is r? Unless the imaginary parts are several orders of magnitude less than the real parts, they may have significant information. Consider using the abs function instead of real if you want to retain that information.) If you’re going to track the relative positions, velocities, altitudes and bearings of all aircraft in the sector, you will probably need to run your loop for each pair of aircraft every minute. (For 12 aircraft, that’s 66 iterations; for 24 aircraft, 276.) You may want to consult the literature to see if someone has devised a much more efficient way of doing this. I’m not aware of one, but I never did such a study.

BTW, commenting your code to remind yourself what it is doing (line-by-line or by groups of lines) is good practice. Comments make it easier for others to read your code — and for you to remember what you did when you come back to it after time away from it.

3. Separation of 500 feet would definitely get my attention!

‘November ###, traffic twelve o’clock, opposite direction, type and altitude unknown’ always provided a little excitement as well.

---------------------------------------------------------------------------

FWIW, I’m an Instrument-Rated Private Pilot (Board Certified Internal Medicine Physician and M.S. Biomedical Engineer), and once designed a study to measure stress in Air Traffic Controllers at the Albuquerque ARTCC. (We had the FAA and NATCA on board and were approved by the NHLBI — highest score — but not funded — too expensive.) Other research I was intimately involved in used EEG to determine what task a person was performing — a step above the cockpit voice recorder! Got to fly an F-16 sim as a research subject in a colleague’s study.

[Caitlin]

Wow - you have a lot of similar interests and experiences! i am also an IFR pilot - and what we are essentially trying to develop here is a function, safety critical task for ATC to use to help preemptively solve critical separation issues.

the k1-k4 equations are essentially functions to combine all the geometry and trigonometry together to find the distance between two aircrafts, subject to the possibility of this "psi" - which psi is if either aircraft were to make an error; pis is the degree of the turn that is made in error. So what it all around calculates is essentially....how many minutes exist until, if either aircraft were to make one of these errors (hence the array), they'll be only 500 ft away from each other?

Does that makes sense? So basically i'm saying...under normal circumstances, they're appropriately separated and its fine. but situations DO occur where if one plane misses its turn, they will almost immediately lose separation with a neighboring plane. so that's what i'm calculating...the time until situations like that COULD occur, so controllers can preemptively solve those issues before they become a VERY serious safety issue.

OK - so that being said. yes, you are correct, i need to evaluate PAIRS of aircrafts. and i thought about it this weekend, and yes..the for loops above aren't really want i want to do..i need to evaluate the pairs.

the for loop is correct for the error vector, the psi vector...that's what i want. for each pair, i need to know the result for each error possible. but then i need to evaluate the PAIRS. so aircraft 1 and 2, 1 and 3, 1 and 4...etc. for each minute. sadly, there is not an efficient way to do this....it's a limitation of my work, these kinds of non-linear functions are extremely complex to optimize.

do you have any ideas on how to evaluate pairs like this? slow and inefficient is OK - i'll have to harass the optimization department at my university to get it any faster..and that's down the line!

i changed all my sin/cos to sind and cosd - that's a huge help, thanks. no one ever taught me that in undergrad (oh so many years ago!)

thanks for your help and collab - and i'd love to read any papers or work you did on your ATC studies and your EEG work for task performance - i do task performance, cockpit design, and function allocation for pilots also! my email is [e-mail] ;-)

thanks@

[Star Strider]

My pleasure!

My ATC stuff never made it to print because we were never able to do that study. (It was also about 25 years ago, so likely obsolete by now.) The EEG stuff is being continued not the group I was with, but by Pfurtscheller and his collaborators. (They and we presented similar papers on essentially the same topics at the same conference in 1992. We published in 1995. That paper is not online and the technology has advanced considerably since.) I would do a PubMed search on "Pfurtscheller[au]" to get those papers. (Shortly before I joined the group I was with, that group had completed the human factors studies and finalized the cockpit design for the C17, probably the first time formal human factors research was ever used in cockpit design. I only read the studies and saw the results.)

I was thinking about your doing multiple simultaneous calculations on multiple aircraft, and the only solution I could come up with is to calculate them as matrices. I’m not quite sure how you would do that with your functions because they appear complicated, and you would likely end up multiplying and adding matrices. It would require more memory and possibly multiple matrices, but you could probably do your calculations much faster as the result. MATLAB has sophisticated matrix functions, and many functions (including the trig functions like cosd) can take matrix arguments, so this would probably be the way I would go. Vectorizing some statements would speed your calculations as well. You can then scan the resulting matrices relatively quickly. Your theta, v, x and y variables would lend themselves to be matrices, to calculate k1 ... k4 and possibly k5 ... k7, you could also use matrices. The problem arises in calculating p efficiently from the matrix entries, so you might not be able to completely eliminate a loop or two. Since you have optimization and MATLAB experts at Purdue, I would consult with them to optimize your code. (That’s not efficiently done online, unfortunately.)

The degree trig functions are relatively new — within the last three or four releases if I remember correctly — so they may not have existed when you initially learned MATLAB. They’re more accurate than converting degrees to radians with pi/180.

Since this is an area of mutual interest and I understand your research, I’ll look around the Internet, and if I see anything that looks promising, I’ll send it along. I’ll help as I can, otherwise. I became interested in human factors science as the result of my getting my IFR rating. Workload saturation and getting vertigo at 12,000 feet in the clouds over the mountains while switching NAV frequencies in an iced airplane focused my interest. (A quick cross-check and an almost nausea-inducing but otherwise uneventful recovery followed.)

[Caitlin]

i thought long and hard yesterday - amongst the distraction of the cat dragging a live baby bunny into the house - and here's what i've come up with.

i'm making one huge matrix, per minute. i've chosen to do that because... the data comes in packets like that anyways. it would, shockingly, be more effort to put it all together than to break it up.

OK so each minute of data has one matrix that is n rows (n=# of aircrafts) and 4 columns; v1 (speed), heading (theta), x, y (positions in space). and then ....

datamatrix=[674 0 6049 46489;760 180 90735 120980;680 0 60049 58500];

psivct=[-90, -45, 0, 45, 90];%degrees degrees of error (should be array)

psidot=3; %degrees/sec rate of error change

d=500; %distance they should be apart for a problem, in feet

for c1 = 1:size(psivct,2)

    psi1=psivct(c1); %for loop to run through array of blunder errors

for i=[1:3]; %to cycle through i=# aircrafts, here #=3

    n=[1,2]; %n=adder to evaluate next aircraft

while n<=i; %i never want it to evalute pairs 1,1 or 2,1...so the # of the first aircraft should always be great than the second.

k1=datamatrix(i,1).*cosd(datamatrix(i,2)+psi1)-datamatrix(i,1)*cosd(datamatrix((i+n),2));

k2=((datamatrix(i,1)/psidot).*(sind(datamatrix(i,2)+psi1)-sind(datamatrix(i,2))))-cosd(datamatrix(i,2)+psi1)*datamatrix(i,1).*(psi1./psidot)+datamatrix(i,3)-datamatrix((i+n),3);

k3=datamatrix(i,1).*sind(datamatrix(i,2)+psi1)-datamatrix((i+n),1)*sind(datamatrix((i+n),2));

k4=(datamatrix(i,1)/psidot)*(cosd(datamatrix(i,2))-cosd(datamatrix(i,2)+psi1))-sind(datamatrix(i,2)+psi1)*datamatrix(i,1).*(psi1./psidot)+datamatrix(i,4)-datamatrix(i+n,4);

k5=k1.^2+k3.^2; %"a" in the polynomial

k6=2*(k1.*k2+k3.*k4); %"b" in the polynomial

k7=k3.^2+k4.^2-500; %"c" in the polynomial

p = [k5 k6 k7]; % creating the polynomial

r = roots(p); %solving for roots

rr(:,c1)=real(r)/60 %getting real answers

end; end; end;

OK - so i'm having issues with my conditions. i'm not sure i'm writing them correctly. essentially, in my formulas, the column that i'm calling for is static, it will never change, so just the row will. so i'm trying to create conditions to say....run the formulas for i=1 to # of aircrafts i have (in this case, 3) and then when i need to call on the second aircraft in the pair, i will add the "+n" to it...so it will evaluate 1 and 1+1=2, 1 and 1+2=3, etc. Do you understand what I mean?

What i'm trying to avoid in unnecessary calculations, so i want to avoid it evaluating 1 and 1, and also evaluating 1 and 2 and then 2 and 1. So that's why i want to always make sure that i > i+n ...does that make sense? I"m not sure i'm coding it right, i'm getting a lot of ...matrix coding errors and things aren't the right dimensions..etc. Let me know what you think...

I'm not much of a programmer - more of a theorist so this stuff is kind of like...oy! and i learned matlab literally in 2005...a LOT has changed in programming since then so i feel totally behind it. I can express it mathematically on paper, but in matlab coding...i get lost :-)

Thanks in advance - I like the dialogue going on...and YIKES - icing is the scariest, and add that in with vertigo, oh man. i've been fairly lucky so far, i've had some scary IFR landings with literally almost no visibility, i lost a landing gear once and lived to tell the tale, but other than that, i've eeked by safely!!!!!

[Star Strider]

One matrix per minute of observations for each variable is likely the best design in your situation. You will likely have to carry over (perhaps correcting) groundspeed and course from one minute to the next.

The easiest way to use loops to avoid duplicates is to have an inner loop counter begin with an increment of the outer loop counter:

for k1 = 1:Nr_AC-1
    for k2 = k1+1:Nr_AC
        AC(k1,k2) = ... % Data and computations for AC#k1 and AC#k2
    end
end

with AC indicating aircraft pairs and Nr_AC being the number of aircraft under consideration. This way, you’ll never end up comparing an airplane with itself. (You’ll end up with a upper-triangular matrix with NaN values on the main diagonal in each variable that uses one, so there may be ways to use that to your computational advantage.) The obvious problem here is that Nr_AC is going to change unpredictably from one minute to the next, as aircraft enter and leave the sector. Probably the easiest (but likely not the most efficient) way to handle this problem is to create an aircraft matrix that is larger than any number you would likely get (perhaps 30), and padding the unused rows and columns with NaN values. You will have to test and eliminate the NaN values at each iteration, but I see no other way to do what you want. You will also have to figure out a way to deal with the aircraft in a specific row or column disappearing from one minute to the next, and figure out where and how to add new aircraft to your matrix. You’ve likely already considered and designed these into your calculations, but are simply items I’ve thought of in the interim.

I’ve been working with MATLAB on and off since 1993, when it was distributed on a stack of 3½" floppies. It’s significantly improved (and gotten easier to work with) over the years, but has changed. (I started programming in college, with FORTRAN II on punchcards.)

Your cat was likely doing it best to impress you and get in your good graces!

I had to drop the gear manually once, but my lowest ceiling was only about 500 feet AGL. A 12-kt gusting crosswind made that landing memorable, though.

[Caitlin]

turns out - i found a way to run through iterations just using a call to a matrix!

c=combk(1:Nr_AC;2)

will give me a matrix c of all possible combinations of aircrafts 1-Nr_AC, and then from that...i can pull into my equations which pairs to reference and then pull the actual data from the data matrix.

done and done!

and as far as the # of aircrafts changing as time goes on in the sector...luckily, this is an experimental data with many limitations and assumptions, so holding the # static at the largest possible number is fine, covering a "worst case scenario."

whew :-)

[Star Strider]

Great!

I’d have suggested combnk but I thought you had to go through all aircraft, not just selected ones. If you’re sampling a representative subset of your aircraft, that will work.

I’m glad you were able to restrict your sector traffic to make your problem more tractable. That problem concerned me, because it’s not obvious.

I don’t know if you intend to publish your research (not all theses and dissertations are published in referred journals), but if you do, I’d appreciate a PDF. My contact information is in my File Exchange profile.

Have fun, stay safe, and avoid occluded fronts!