0
$\begingroup$

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

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

enter image description here

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

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

enter image description here

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

enter image description here

Finally get h3 from the third equation

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

enter image description here

Do the same with the second solution of sol1

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.