# Contour Plot of a set of equations? [closed]

For the set of equations the set of equations, w versus A has to be plotted, I have written the code but can't seem to figure out where I'm wrong, any help would be appreciated. Thanks.

Clear[c, m, e, c1, cp, M, F, Fp, A1, B1, A1p, B1p, bet, alp, betp, \
alp, w, wpp]
m = 0.08
e = 0.4
Mum=0.08
Rm=1
c1 = 0.02
w =.
A =.
wpp =.
k1 = 1000
k2 = 700
lep = 3*(k1/(M + m))^0.5
M = 1
c = c1/(M + m)
const1 - 4*funct2
alp = k1/(M + m)
bet = k2/(M + m)
F = m*e*w^2/(M + m)
wpp = w/alp^0.5
cp = c/(alp)^0.5
betp = bet*A^2/alp
Fp = F/(alp*A)
A1p = A1/A
B1p = B1/A
eqn = {-cp*wpp/3 + 3*betp*B1p/4 == 0,
((1 - wpp^2)*B1p + (3*betp*B1p/4)*(2 + A1p^2 + B1p^2) -
cp*wpp*A1p)^2 + ((1 - wpp^2)*A1p + (betp/4)*(1 + 6*A1p + 3*A1p^3 +
3*A1p*B1p^2)+cp*wpp*B1p)^2 == Fp^2,
(1 - wpp^2/9) + (3*betp/4)*(1 + A1p + 2*A1p^2 + 2*B1p^2) ==0,
Mum*Vs/Rm - Mum^2*w/Rm == c*(A^2*w^2/18 + Ah^2*w^2/2),
((alp - w^2)*Ah + 3/4*bet*Ah^3)^2 + c^2*Ah^2*w^2 == (F*w^2)^2
};
ContourPlot[eqn, {Vs,0,15},{w}]

• B1 can be found out in terms of A from -cpwpp/3 + 3*betpB1p/4 == 0 and A1 can be found out from ((1 - wpp^2)*B1p + (3*betpB1p/4)*(2 + A1p^2 + B1p^2) - cpwppA1p)^2 + ((1 - wpp^2) A1p + (betp/4)*(1 + 6*A1p + 3*A1p^3 + 3*A1pB1p^2) + cpwpp*B1p)^2 == Fp^2 in terms of A, basically reducing the whole set into a relation between A and w and plotting for {w,90,95} particularly where the subharmonic amplitude exists. – Pavan Kumar Apr 12 '17 at 16:16
• edit the question with new info. Its too hard to read code in comments. – george2079 Apr 12 '17 at 18:34
• Edited, Please have a look. – Pavan Kumar Apr 12 '17 at 20:35
• You have significantly changed the equations in the question. You now have eight unknowns and five equations. What do you want plotted against what? – bbgodfrey Apr 12 '17 at 22:15
• @bbgodfrey The first three equations in eqn give a relationship between A and w,The fourth is a relation between Vs,A and Ah, the final one is a relation between Ah and w, the plot is for Vs and w. for {Vs,0,15} – Pavan Kumar Apr 12 '17 at 22:34

The system of equations eqn (less the second equation) is equivalent to a 24th order polynomial, so eliminating variables {A, A1, Ah, B1} symbolically does not seem practical. Instead, eliminate those variables numerically, proceeding as follows. First, Rationalize all decimal numbers (perhaps unnecessary here but usually a good practice). Then define a new variable f as the Real numerical mismatch between the left and right sides of the second equation. It is this quantity that will be used in ContourPlot.

eqnr = Rationalize@eqn;
eqnf = Join[eqnr, {f ∈ Reals}]; eqnf[[2]] = f == Subtract @@ eqnr[[2]]


and actually evaluate f by means of

ff[Vs0_?NumericQ, w0_?NumericQ] :=
Min[f /. NSolve[eqnf /. {Vs -> Vs0, w -> w0}, {A, A1, Ah, B1, f}]]


which solves for f numerically, and returns the smallest value, if any. A typical plot of ff for constant Vs is

Plot[ff[1., x], {x, 9, 12}, PlotRange -> {Automatic, {-.02, .2}}, PlotPoints -> 10]


It shows three roots, the second of which can be traced as a function of Vs in a few minutes.

opts = {AccuracyGoal -> 4, PrecisionGoal -> Infinity}
start = {x /. Quiet@FindRoot[ff[1, x], {x, 9.7, 9, 10}, Evaluate@opts],
x /. Quiet@FindRoot[ff[1.1, x], {x, 9.9, 9, 10}, Evaluate@opts]};
Do[tem = 2 start[[i - 1]] - start[[i - 2]];
new = x /. Quiet@FindRoot[ff[1. + (i - 1)/10, x], {x, tem, tem - 0.3, tem + .3},
Evaluate@opts]; If[NumericQ[new], AppendTo[start, new], Break[]], {i, 3, 141}]
p1 = ListLinePlot[start, DataRange -> {1, 1 + (Length@start - 1)/10}]


Finally, the requested ContourPlot can be computed, although slowly, by

Table[{y, x, ff[y, x]}, {y, 1, 15, .02}, {x, 9, 14, .02}];
p2 = ListContourPlot[Flatten[%, 1], Contours -> {10^-3}, ContourShading -> None,
AspectRatio -> 1/GoldenRatio]


The "zero" contour has been chosen as slightly larger than zero, because the value of ff never quite reach zero, because of the finite precision of the original equations in the question as well as the finite resolution of ListContourPlot. Even with this chose of contour value, the curves are broken in places.

here is a start, make all your floats rationals (eg, m=8/100) and note the first equation is linear in B1, so eliminate to obtain two equations in A,A1,w:

 eqn = Simplify[(eqn /. First@Solve[eqn[[1]], B1])[[2 ;;]]]


ContourPlot3D can handle this:

ContourPlot3D[Evaluate[eqn], {w, 0, 150}, {A, -2, 2}, {A1, -4, 4},
PlotPoints -> 10]


the intersection of those surfaces is your solution ( so at least we know there is a solution.

Now you can play the same game and eliminate A1 , the last equation is quadratic so you get two solutions:

  {v1, v2} = eqn[[1]] /. Solve[eqn[[2]], A1];


unfortunately this doesn't work..probably because the function produces many complex values near the zero solution

  ContourPlot[v1, {w, 100, 120}, {A, -2, 2}]


looking at the real part:

g = Evaluate[v1[[1]] - v2[[2]]];
ContourPlot[Re[g] == 0, {w, 0, 150}, {A, -2, 2}, Contours -> {0}]