# Three nonlinear equations with three unknowns and plot one of the unknowns with a parameter?

There are three nonlinear equations of unknowns h1, h2 and h3, and I need to plot the unknowns with varying $\sigma_1$.

1.*10^-6 ==0.0001 h1^2 + 100 h1 ((51 h2^2)/1400 - Subscript[σ, 1])^2

1.*10^-6 ==4.*10^-6 h2^2 + 4 h2^2 (0.00104167 + (6 h1^2)/7 - (33 h2^2)/350 -
Subscript[σ, 1])^2

1.*10^-6 ==1.6*10^-7 h3^2 + 0.16 h3^2 (0.001 + 0.00462963 h1^2 + 0.00204167
h2^2 - 0.0125 h3^2)^2

• Have you tried anything? – Kuba Feb 16 '18 at 15:40
• Try FindRoot. – Henrik Schumacher Feb 16 '18 at 16:10
• What is the range of values for $\sigma_1$? – Carl Woll Feb 16 '18 at 16:15
• the range of sigma 1 is from -0.008 to 0.008 – Mohamed Khalid Feb 16 '18 at 18:45
• @Henrik Schumacher how can i use FindRoot – Mohamed Khalid Feb 16 '18 at 20:01

For problems like this, I like to recast it as an ODE where the parameter is the independent variable. Here is how I would do this for your question.

Equations

First, make the dependence of the dependent variables explicit:

eqns = {
1.*10^-6==0.0001 h1^2+100 h1 ((51 h2^2)/1400-Subscript[σ,1])^2,
1.*10^-6==4.*10^-6 h2^2+4 h2^2 (0.00104167+(6 h1^2)/7-(33 h2^2)/350-Subscript[σ,1])^2,
1.*10^-6==1.6*10^-7 h3^2+0.16 h3^2 (0.001+0.00462963 h1^2+0.00204167 h2^2-0.0125 h3^2)^2
};

peqns = eqns /. {h1->h1[s], h2->h2[s], h3->h3[s], Subscript[σ,1]->s};
peqns //TeXForm


$\left\{\text{1.$\grave{ }$*${}^{\wedge}$-6}=100 \operatorname{h1}(s) \left(\frac{51 \operatorname{h2}(s)^2}{1400}-s\right)^2+0.0001 \operatorname{h1}(s)^2,\text{1.$\grave{ }$*${}^{\wedge}$-6}=4 \operatorname{h2}(s)^2 \left(\frac{6 \operatorname{h1}(s)^2}{7}-\frac{33 \operatorname{h2}(s)^2}{350}-s+0.00104167\right)^2+\text{4.$\grave{ }$*${}^{\wedge}$-6} \operatorname{h2}(s)^2,\text{1.$\grave{ }$*${}^{\wedge}$-6}=0.16 \operatorname{h3}(s)^2 \left(0.00462963 \operatorname{h1}(s)^2+0.00204167 \operatorname{h2}(s)^2-0.0125 \operatorname{h3}(s)^2+0.001\right)^2+\text{1.6$\grave{ }$*${}^{\wedge}$-7} \operatorname{h3}(s)^2\right\}$

Initial equations

Next, to use NDSolveValue, we need an initial condition, which we find using NSolve:

p0 = First @ NSolve[(And@@peqns /. s->1) && h2[1]>0 && h3[1]>0, Reals]


{h1[1] -> 1.*10^-8, h2[1] -> 0.000500521, h3[1] -> 0.624558}

There are many possible solutions, so I chose the one with positive numbers.

In versions before M10, you may need to use FindRoot instead of NSolve:

FindRoot[peqns/.s->1, {{h1[1],1}, {h2[1],1}, {h3[1],1}}]


{h1[1] -> 1.*10^-8, h2[1] -> 0.000500521, h3[1] -> 0.624558}

NDSolveValue

Now, we can use NDSolveValue:

sol = NDSolveValue[
Join[D[peqns, s], Equal @@@ p0],
{h1, h2, h3},
{s, 0, 5}
];


NDSolveValue::ndsz: At s == 0.001652735639801279, step size is effectively zero; singularity or stiff system suspected.

Visualization

GraphicsRow[{
Plot[sol[[1]][s], {s, 0.002, 5}],
Plot[sol[[2]][s], {s, 0.002, 5}],
Plot[sol[[3]][s], {s, 0.002, 5}]
}]


• this code give me error from po= First..... – Mohamed Khalid Feb 16 '18 at 19:17
• is there another way to solve this problem – Mohamed Khalid Feb 16 '18 at 19:23
• @MohamedKhalid You probably have lingering definitions. You could either quit the kernel and retry, or do ClearAll["Global*"] and retry. – Carl Woll Feb 16 '18 at 19:23
• i tried quit kernel and clear all and i have the same error – Mohamed Khalid Feb 16 '18 at 19:37
• I do two orders but i have the same error – Mohamed Khalid Feb 16 '18 at 19:56

Here a solution by solving one equation after the other.

First rationalize equations

eqs = Rationalize[eqns, 0] /. Subscript[\[Sigma], 1] -> s1;


Solve the first equation for h1, since it does not depend on h3

sol1 = Solve[eqs[[1]], h1]

(* {{h1 -> (1/
9800)(-6502500 h2^4 + 357000000 h2^2 s1 - 4900000000 s1^2 - Sqrt[
960400 + (-6502500 h2^4 + 357000000 h2^2 s1 -
4900000000 s1^2)^2])}, {h1 -> (1/
9800)(-6502500 h2^4 + 357000000 h2^2 s1 - 4900000000 s1^2 + Sqrt[
960400 + (-6502500 h2^4 + 357000000 h2^2 s1 -
4900000000 s1^2)^2])}}   *)


You get two solutions. Now insert one of that solutions into second equation and use FindRoot to get the dependance of h2 on s1.

fh2[s_] :=
First@FindRoot[
eqs[[2]] /. sol1[[1]] /. s1 -> Rationalize[s, 0], {h2, .01},
WorkingPrecision -> 200, AccuracyGoal -> 20, PrecisionGoal -> 20]

LogPlot[h2 /. fh2[s], {s, -.008, .008}]


Get h1 now from the first equation

fh1[s_] := First@FindRoot[eqs[[1]] /. fh2[s] /. s1 -> s, {h1, .05}]

LogPlot[h1 /. fh1[s], {s, -.008, .008}]


Finally get h3 from the third equation

Plot[h3 /.
First@FindRoot[
eqs[[3]] /. fh2[s] /. fh1[s], {h3, .5}], {s, -.008, .008},
PlotRange -> All]
`

Do the same with the second solution of sol1