# A system of 2 non-linear equations with 2 unknowns

1. Nov 4, 2013

### pm1366

1. The problem statement, all variables and given/known data

i'm trying to find a solution for these two equations, p & q are variables and c is known constant (it's given randomly) :

\begin{align} F_1=&\left(p^2q^2(p^2-q^2)+c^2(p^4+q^4-6p^2q^2)-3c^4(p^2-q^2)-4c^6\right)\sin(p)\sinh(q)+\\&2pq\left(-p^2q^2+3c^2(p^2-q^2)-7c^4\right)\cos(p)\cosh(q)-\\&pq\left(p^4+q^4+2c^2(q^2-p^2)+2c^4\right)\cos(2c) \\ F_2=&2cp\left(3p^2q^2-q^4-c^2(p^2+q^2)+4c^4\right)\cos(p)\sinh(q)+\\&2cq\left(p^4-3p^2q^2-c^2(p^2+q^2)-4c^4\right)\sin(p)\cosh(q)+\\&pq(p^2+q^2)(q^2-p^2+2c^2)\sin(2c) \end{align}

(F1 is taken from real part of an equation and F2 is taken from its imaginary part and both are equal to zero)

i want to solve F1 & F2 to find p and q for a given c

p and q are real variables (they were taken from a complex unknown previously p+qi, i mean we know that p and q must be Real)

2. Relevant equations
all above

3. The attempt at a solution
i tried with both maple and matlab with solve command but it takes a long time running for them and finally there's no answer & my PC gets slow !

any procedure or solution that can help or maybe a way to solve with maple or matlab is appreciated .

2. Nov 4, 2013

### UltrafastPED

If you know the range of valid values for c you can solve the equations numerically:
do one value of c at a time, and plot the p,q that is obtained. Start with large step sizes!

3. Nov 4, 2013

### Ray Vickson

Your problem is unclear: are $F_1$ and $F_2$ given numerical values, or are the just shorthand names for the formulas on the right. If they are the latter, then you need to specify what numerical values they should have.

In Maple you can try an "fsolve" command, and perhaps specify ranges on p and q if you can determine them somehow. In other words, if F1 and F2 are given values, and f1, f2 are the functions of p and q you wrote above, you can try the command

fsolve({f1=F1,f2=F2},{p,q}).

You might need to do some 3d plots first in order to get a feel for the region containing the solution, then you can use the ranges as an option in the fsolve command.

If that does not work you can try instead to use the Optimization package in Maple, to tackle the min-square deviation problem:

with(Optimization);
NLPSolve((f1-F1)^2 + (f2-F2)^2)

Here, if you have inputted numbers for c, F1 and F2, Maple will understand that the only variables are p and q. Again, if you have some rough information about the location of p and q you can incorporate that using the option {p=pmin..pmax, q=qmin..qmax}. The optimization default is minimization, so you need not specify that.

Note that you could use the solution of the minimization problem as input to the equation-solving problem, in order to increase the accuracy of the final solution.

4. Nov 5, 2013

### pm1366

thank you very much Ray , very helpful Notes , i used fsolve and it is ok now ,giving me good values for p&q ,you did a great job , it would be very nice to help me with matlab too in my next comment

Last edited: Nov 5, 2013
5. Nov 5, 2013

### pm1366

what i have is F1=0 & F2=0 , i give some values to c in range of [0,2.5] and fsolve in maple is ok ,this way as you said :

fsolve({F1=0,F2=0},{p,q})

but as i'm not too much familiar to maple , i need the similar case for matlab to be able to write my code easier (matlab is easier for me to understand)
i want to use fsolve in matlab but i get several errors , fsolve in matlab needs an initial point x0 , and i don't understand it well ! it also needs to pass through a function

in a separate m file , i write:

function F=sys(p,q)
global c
F=[((p^2)*(q^2)*((q^2)-(p^2))+(c^2)*((p^4)+(q^4)-6*(p^2)*(q^2))-3*(c^4)*((p^2)-(q^2))-4*(c^6))*sin(p)*sinh(q)+2*p*q*((-(p^2)*(q^2))+3*(c^2)*((p^2)-(q^2))-7*(c^4))*cos(p)*cosh(q)-(p*q)*((p^4)+(q^4)+2*(c^2)*((q^2)-(p^2))+2*(c^4))*cos(2*c); 2*c*p*(3*(p^2)*(q^2)-(q^4)-(c^2)*((p^2)+(q^2))+(4*(c^4)))*cos(p)*sinh(q)+2*c*q*((p^4)-3*(q^2)*(p^2)-(c^2)*((p^2)+(q^2))-4*(c^4))*sin(p)*cosh(q)+(p*q)*(((p^2)+(q^2))*((q^2)-(p^2)+2*(c^2)))*sin(2*c)]
end

and in another m file :

global c
c=1; %for example
fsolve(@sys)

But i get this error:
??? Error using ==> fsolve at 138
The input to FSOLVE should be either a structure with valid fields or consist of at least two
arguments.

Error in ==> test2 at 6
fsolve(@sys)

6. Nov 5, 2013

### Ray Vickson

I have not used Matlab, and do not have access to it.

7. Nov 5, 2013

### pm1366

no problem Mr Vickson; there is a contradiction here ! when i change the range of p and q i get different results :
for example ,if i set :
> S := fsolve({F1 = 0, F2 = 0}, {p = 0 .. 20, q = 0 .. 20});

{p = 15.03019936, q = 0.}

and if :
> S := fsolve({F1 = 0, F2 = 0}, {p = 0 .. 10, q = 0 .. 10});

{p = 8.031885010, q = 0.}

and if i set no ranges , this way:
> S := fsolve({F1 = 0, F2 = 0}, {p, q});

{p = 7.539612454, q = -5.355194539}

i can't find the correct range , and i don't know how to 3d plot in maple , what did you mean when you said 3dplot ? 3d plot of what ? can you make it clear for me? i'm newbie in maple . thanks

Last edited: Nov 5, 2013
8. Nov 5, 2013

### UltrafastPED

3D plot of p,q,c.

9. Nov 5, 2013

### pm1366

thanks but i am seeking for p,q , how can i plot them when don't have them in hand ?! i want to know a possible range for solution of p & q and then apply the range to fsolve command ...

10. Nov 5, 2013

### UltrafastPED

I cannot answer your Matlab questions; I am just suggesting generally applicable mathematical ideas.

For example, for a given value of c there should be a set of p,q values, each corresponding to one complex number z. That is, for each c you should obtain a curve on the complex plane.

As you solve for additional values of c you will generate additional curves; this family of curves is parameterized by the value of c.

I don't think Matlab likes to work with complex numbers ... you might do better to use Mathematica or Maple.

11. Nov 5, 2013

### pm1366

^ i'm now considering only maple , what i want to know is a suitable range for fsolve and just for 1 value of c , say c=1 for example ,
in maple for a specific value of c :

> S := fsolve({F1 = 0, F2 = 0}, {p = 0 .. 20, q = 0 .. 20});

{p = 15.03019936, q = 0.}

and if :
> S := fsolve({F1 = 0, F2 = 0}, {p = 0 .. 10, q = 0 .. 10});

{p = 8.031885010, q = 0.}

and if i set no ranges , this way:
> S := fsolve({F1 = 0, F2 = 0}, {p, q});

{p = 7.539612454, q = -5.355194539}

as you see with every range for p &q , i get different results , now i need a range that make sense , a good range for p & q which give me the best result ...

12. Nov 5, 2013

### Ray Vickson

I meant a plot of the surface z = f1(p,q) and z = f2(p,q), but I now think that implicitplots of f1(p,q) = 0 and f2(p,q)=0 are more appropriate. For example, when c = 0 you can put f10:=subs(c=0,f1) and f20:=subs(c=0,f2), then use
with(plots);
[animate, animate3d, animatecurve, arrow, changecoords, complexplot,

complexplot3d, conformal, conformal3d, contourplot,

contourplot3d, coordplot, coordplot3d, densityplot, display,

implicitplot, implicitplot3d, inequal, interactive,

interactiveparams, intersectplot, listcontplot,

listcontplot3d, listdensityplot, listplot, listplot3d,

loglogplot, logplot, matrixplot, multiple, odeplot, pareto,

plotcompare, pointplot, pointplot3d, polarplot, polygonplot,

polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus,

semilogplot, setcolors, setoptions, setoptions3d, spacecurve,

sparsematrixplot, surfdata, textplot, textplot3d, tubeplot] <--- a list of all the plot forms

implicitplot(f10=0,p=-20..20,q=-10..10);
The results are not very good (they are improvable, but that does not matter much here); you will see that the points (p,q) satisfying f10(p,q) = 0 consist of the p and q axes, plus numerous disjoint vertical sausage-shaped curves above and below the p axis. The points satisfying f20(p,q)=0 are the whole (p,q) plane, since f20(p,q) = 0 for all p and q. Therefore, there are infinitely many solutions when c = 0.

For c = 0.1 you can do f11:=subs(c=0.1,f1);f21:=subs(c=0.1,f2);
G1:=implicitplot(f11=0,p=-20..20,q=-10..10,color=red): <--- note colon, not semi-colon
G2:=implicitplot(f21=0,p=-20..20,q=-10..10,color=blue):
display([G1,G2]); <--- this overlays the two graphs together
While the implicitplot results are, again, not very good, they are good enough to show what is happening: there are many, many solutions to the nonlinear system, corresponding to all the points where the red and blue curves cross.

Therefore, it is not at all surprising that you do not get a unique solution; there must be some other information you need to supply in order to get a unique solution, or at least a smaller number of possible solutions.

I am assuming your formulas for F1 and F2 are correct; I just displayed them as tex formulas, then copied them to a Maple worksheet and edited them from tex form to Maple
form.

Last edited: Nov 5, 2013
13. Nov 6, 2013

### pm1366

thanks Mr Vickson ,
yes i saw those plots , it seems that this system has a lot solutions , [ i made a simple mistake in F1 , instead of (p^2 - q^2) we must have (q^2 - p^2) in first term , although this doesn't affect the whole problem and plots ]

so the correct F1 & F2 are:

\begin{align} F_1=&\left(p^2q^2(q^2-p^2)+c^2(p^4+q^4-6p^2q^2)-3c^4(p^2-q^2)-4c^6\right)\sin(p)\sinh(q)+\\&2pq\left(-p^2q^2+3c^2(p^2-q^2)-7c^4\right)\cos(p)\cosh(q)-\\&pq\left(p^4+q^4+2c^2(q^2-p^2)+2c^4\right)\cos(2c) \\ F_2=&2cp\left(3p^2q^2-q^4-c^2(p^2+q^2)+4c^4\right)\cos(p)\sinh(q)+\\&2cq\left(p^4-3p^2q^2-c^2(p^2+q^2)-4c^4\right)\sin(p)\cosh(q)+\\&pq(p^2+q^2)(q^2-p^2+2c^2)\sin(2c) \end{align}

to tell you the whole story of the problem , i must say that after obtaining a solution for p & q , we have these 3 extra equations:

$$F_3= u^2-p^2+q^2-2*c^2$$

$$F_4= ω^2-(p^2-c^2)*(q^2-c^2)$$

$$F_5= \beta ^{0.5}*u*ω-c*(p^2+q^2)$$

then we can solve these 3 equations with 3 unknowns to find ω & β & u
also we know that β≠0 & c≠0

the whole computational procedure is as follows:

1- select a Real value for c
2- finding values of p & q satisfying F1=0 and F2=0 in a system
3- determine u and ω and β from F3 ,F4 and F5

the only thing that i know is that the resulting β must be between 0 and 1 , i mean the range (0,1)
i mean we must first select the values of c in a way to cover final β in the range (0,1)
we must repeat the above procedure to cover all values of beta, then we must plot u-β and ω-β figures

i hope you can help me , all your helps are appreciated

Last edited: Nov 6, 2013
14. Nov 6, 2013

### Bill Simpson

IF I haven't made any mistakes then

Using your revised F1 and F2, arbitrarily picking c=1/2 as a first example, letting Mathematica minimize F1^2+F2^2 by randomly picking points and going downhill from each point almost instantly gives one solution of p=1, q=0. Choosing c=1/4 gives p=1/2, q=0. Choosing c=3/4 gives p=0, q=0. etc.

Now I assume I need to know the constants F3, F4 and F5 so that I can solve those three equations in 3 unknowns. Or perhaps those are also zero.

If you have a range of "interesting" c values that would be helpful to know. Or if you really trying to find all acceptable c values with the only condition being 0<beta<1 that would be helpful to know.

If I assume for the moment F3=F4=F5=0, p=1,q=0,c=1/2 then ω^2= -3/16. Are complex solutions acceptable? Or are there more conditions? Otherwise this solution is
u= -Sqrt[3/2], ω= I Sqrt[3]/4, β= -8/9
or
u= Sqrt[3/2], ω= -I Sqrt[3]/4, β= -8/9
and that doesn't satisfy 0<β<1 so something has to change.

Assuming F3=F4=F5=0 it appears that solutions for small positive rational values of c and resulting non-negative β are all p=q=0 and that gives β=0. So something else needs to reject p=q=0.

Last edited: Nov 6, 2013
15. Nov 7, 2013

### pm1366

thanks Mr Simpson for your help , but i get different results for solution of the system with Maple:

Even if i give positive range to p & q before solving , this way:
S := fsolve({F1 = 0, F2 = 0}, {p = 0 .. 10, q = 0 .. 10})

and supposing c=1/2 , i get {p = 7.991211001, q = 0.}

and also with the same c, if i give a range between -10 , 10 , this way:
S := fsolve({F1 = 0, F2 = 0}, {p = -10 .. 10, q = -10 .. 10})

then the result is {p = 5.255360827, q = -2.697303464}

which are toally different from yours, unfotunately i know nothing about Mathematica.

with my maple code (with the aid of fsolve) and comparing to a presolved Paper which i'm working on , for beta<0.3 i get good results (similar results with paper) but for β>0.3 it has different values but with a similar shape in the whole figure,

Friends, is there any other Way to solve F1=0 & F2=0 except fsolve command ???

16. Nov 7, 2013

### Bill Simpson

Probably any set of equations involving Sin and Cos are going to have multiple solutions, perhaps even an infinite number of solutions, especially when I used a random search process to stumble onto a solution. It might be a miracle if you happened to get the same solution that I did.

Unfortunately the only copy of Maple I have is extremely old and I don't think it would run on my current system. But a quick Google search for Maple minimize claims that Maple can find a local minimum using the minimize function. Perhaps you can do something similar to what I did, choose ten or ten thousand random starting values for p and q, or walk across an evenly spaced grid of 100x100 p,q values, and use the Maple minimize function on f1^2+f2^2 to look for solutions. Collect all those solutions together, eliminate duplicates and see if you observe any pattern in the results.

You didn't answer whether you then wanted to solve f3=f4=f5=0 or is it something else? Is there any restriction, other than being Real and 0<β<1, on the value of c?

You actually have a published paper that explains how this is all supposed to work... and you haven't told any of us where that paper is??? Why???

Code (Text):

Step through values of c for 1/16<=c<=5/2,
Randomly search for p,q 1/8<p,q<10 to solve f1=f2=0,
Find u,w,b such that f3=f4=f5=0 and 0<b<1.
p        q        c       u        w        b
{{5.30904, 2.64458, 0.56250, 4.67172, 13.6417, 0.096414},
{5.34017, 2.61485, 0.59375, 4.73128, 13.5147, 0.107779},
{5.37466, 2.58262, 0.62500, 4.79566, 13.3768, 0.120008},
{5.41303, 2.54765, 0.65625, 4.86536, 13.2268, 0.133216},
{5.45592, 2.50962, 0.68750, 4.94107, 13.0635, 0.147557},
{5.50420, 2.46815, 0.71875, 5.02371, 12.8851, 0.163247},
{5.69579, 2.31713, 0.81250, 5.32853, 12.2335, 0.222110},
{5.78379, 2.25472, 0.84375, 5.45823, 11.9639, 0.247919},
{5.89225, 2.18352, 0.87500, 5.61089, 11.6570, 0.279047},
{6.03224, 2.10058, 0.90625, 5.79810, 11.3015, 0.318402},
{6.22623, 2.00151, 0.93750, 6.04299, 10.8848, 0.371638},
{6.52188, 1.88413, 0.96875, 6.39234, 10.4225, 0.449032},
{6.94954, 1.78428, 1.00000, 6.86385, 10.1626, 0.544663},
{7.33651, 1.76952, 1.03125, 7.26774, 10.4449, 0.598678},
{7.60547, 1.80913, 1.06250, 7.53844, 11.0271, 0.610213},
{7.79843, 1.86896, 1.09375, 7.72755, 11.7017, 0.605053}}

Last edited: Nov 7, 2013
17. Nov 8, 2013

### pm1366

thank you for your help, actually that is an old paperback in my university library , it doesn't have much explanations, it's a domestic paper...

now that i see your table, your results are the same as mine in maple when p,q≥1/8
exactly the same except that my c is up to max 1.675 to reach b<1

to answer your question , i should say yes , then i need to solve F3=F4=F5=0 to find beta and omega and u. and also there's no other restriction other than being Real and 0<β<1

as you see in your above table , ω is decreasing until b=0.598678 and then it is increasing , in paper this change is at about b=0.3 , that's why i'm looking for a more exact method for this problem.

also another part of this problem (for getting omega-u plot which needs a minimization that i don't have any idea about it !) ,i've asked here: