5
$\begingroup$

The mathematical problem is shown here. I tried to solve writing the 

A0 = 2.0377638272727268`;
A1 = -7.105521894545453`;
A2 = 9.234000147272726`;
A3 = -5.302489919999999`;
A4 = 1.1362478399999998`;
h0 = 45.5;
σM = 0.00592251
λ1 = 1.025;
λ2 = 1.308;

f1[y_, x_] = A1 + 2 A2 y[x] + 3 A3 y[x]^2 + 4 A4 y[x]^3
b[y_, x_] = h0^2/12 (5 A1 + 8 A2 y[x] + 9 A3 y[x]^2 + 8 A4 y[x]^3)/y[x]^6
g[y_, x_] = -(h0^2/12) (A1 + 2 A2 y[x] + 3 A3 y[x]^2 + 4 A4 y[x]^3)/y[x]^5

Sys = First@NDSolve[{
f1[y, x] + b[y, x] y'[x]^2 + g y''[x] == σM,
y[0] == λ1,
y[15] == λ2
}, y, {x, 0, 15}]

When I evaluate this code, Mathematica answers with the following errors:

Power::infy: Infinite expression 1/0.^6 encountered. >> Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

Power::infy: Infinite expression 1/0.^6 encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

Power::infy: Infinite expression 1/0.^7 encountered. >>

General::stop: Further output of Power::infy will be suppressed during this calculation. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

General::stop: Further output of Infinity::indet will be suppressed during this calculation. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`. >>

As suggested here (section: "StartingInitialConditions") I tried to wrote:

Sys = First@NDSolve[{
      f1[y, x] + b[y, x] y'[x]^2 + g y''[x] == σM,
      y[0] == λ1,
      y[15] == λ2}, y, {x, 0, 15}, 
      Method -> {"Shooting","StartingInitialConditions" -> {y[0] == λ1, y'[0] == 0.1}}]

followed by the following error:

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`. >>

How can I overtake this errors? What am I doing something wrong?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
6
$\begingroup$

There's a g in your differential equation, I think you actually mean g[y, x], right? This is the root of your problem. By the way, why do you make y an argument for functions f1, b, g? (It doesn't hurt though…) In fact, I think it will be conciser not to use function definitions in this case.

Then, your "StartingInitialConditions" for Shooting isn't so proper, I found a better one:

A0 = 2.0377638272727268`;
A1 = -7.105521894545453`;
A2 = 9.234000147272726`;
A3 = -5.302489919999999`;
A4 = 1.1362478399999998`;
h0 = 45.5;
σM = 0.00592251;
λ1 = 1.025;
λ2 = 1.308;

f1 = A1 + 2 A2 y[x] + 3 A3 y[x]^2 + 4 A4 y[x]^3;
b = h0^2/12 (5 A1 + 8 A2 y[x] + 9 A3 y[x]^2 + 8 A4 y[x]^3)/y[x]^6;
g = -(h0^2/12) (A1 + 2 A2 y[x] + 3 A3 y[x]^2 + 4 A4 y[x]^3)/y[x]^5;
eqn = f1 + b y'[x]^2 + g y''[x] - σM;

sol = NDSolve[{eqn == 0 , y[0] == λ1, y[15] == λ2}, y, {x, 0, 15}, 
               Method -> {"Shooting", "StartingInitialConditions" -> 
                          {y[15] == λ2, y'[15] == 0}}]
Plot[y[x] /. sol, {x, 0, 15}]
Plot[eqn /. sol, {x, 0, 15}, PlotRange -> All]

enter image description here

enter image description here

Edit

Some warnings will be generated when running the code above, but as shown in the second picture, the error of the solution is quite small, the only defect is, it doesn't match the left BC 囧. (Sorry I didn't notice it yesterday. ) But now I've got a really proper initial condition:

sol = NDSolve[{eqn == 0 , y[0] == λ1, y[15] == λ2}, y, {x, 0, 15}, 
               Method -> {"Shooting", "StartingInitialConditions" -> 
                          {y[15] == λ2, y'[15] == 1/10 + 1/100 + 2/1000}}]
Plot[y[x] /. sol, {x, 0, 15}, PlotRange -> All, 
     Epilog -> {PointSize[Large], Point[{{0, λ1}, {15, λ2}}]}]
Plot[eqn /. sol, {x, 0, 15}, PlotRange -> All]
(y[0] /. sol) - λ1

enter image description here

enter image description here

{4.45117*10^-6}

I still found it by trial and error (with the help of Table).

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks very much: it works very well!! But why does not it work if as starting point in the Shooting option I set y[0]==lambda1? $\endgroup$
    – Petrus
    Commented May 5, 2013 at 14:09
  • $\begingroup$ @Petrus I found the initial condition above by trial and error… I think it's part of the property of shooting method… sorry I can't give a deeper explain, but now I notice that maybe my answer can be better, I may improve it tomorrow. $\endgroup$
    – xzczd ♦
    Commented May 5, 2013 at 16:21
  • $\begingroup$ @xzxzd thanks for your help $\endgroup$
    – Petrus
    Commented May 5, 2013 at 20:53
  • $\begingroup$ @Petrus See my edit. $\endgroup$
    – xzczd ♦
    Commented May 6, 2013 at 4:49
  • $\begingroup$ xzczd:how the use of Table can help you? $\endgroup$
    – Petrus
    Commented May 7, 2013 at 1:42
0
$\begingroup$

Now I tried to change the interval as follows:

A0 = 2.0377638272727268`;
A1 = -7.105521894545453`;
A2 = 9.234000147272726`;
A3 = -5.302489919999999`;
A4 = 1.1362478399999998`;
h0 = 45.5;
σM = 0.00592251;
λ1 = 1.0253896074561006`;
λ2 = 1.3079437258774012`;

f1 = A1 + 2 A2 y[x] + 3 A3 y[x]^2 + 4 A4 y[x]^3;
b = h0^2/12 (5 A1 + 8 A2 y[x] + 9 A3 y[x]^2 + 8 A4 y[x]^3)/y[x]^6;
g = -(h0^2/12) (A1 + 2 A2 y[x] + 3 A3 y[x]^2 + 4 A4 y[x]^3)/y[x]^5;
eqn = f1 + b y'[x]^2 + g y''[x] - σM;

LA = 0;
LB = 15;

sol = NDSolve[{
     eqn == 0,
     y[LA] == λ1,
     y[LB] == λ2
     }, y, {x, 0, 15}, 
     Method -> {"Shooting", "StartingInitialConditions" -> {y[LB] == λ2, y'[LB] == 0}}
     ]

Show[
     Plot[y[x] /. sol, {x, 0, LB}, PlotRange -> {Automatic, {1, 1.32}}, Frame -> True],
     Graphics[{
     Red,
     Point[{{LA, λ1}, {LB, λ2}}]
      }]
    ]

If I set LB different from 15, something strange happens.. Moreover, even if LB = 15, the boundary condition in x = 0 is not satisfied. Where is the problem?

Improve this answer
$\endgroup$
5
  • $\begingroup$ @I think it's better to add this to your question rather than post it as an answer… I recall a excellent answer by someone in this site to solve the problem here, but I'd like to go to bed now, you can look for it by yourself if you don't want to wait :), it's also a question about shooting method. $\endgroup$
    – xzczd ♦
    Commented May 5, 2013 at 16:18
  • $\begingroup$ @xzczd: Now I'll look for this. How written here, the real BC are $$ \lambda (x\to - \infty) = \lambda_1 \qquad \lambda (x\to \infty) = \lambda_2 $$ I'm trying to have a larger domain for this reason $\endgroup$
    – Petrus
    Commented May 5, 2013 at 20:48
  • $\begingroup$ I found that question: mathematica.stackexchange.com/questions/18301/… , though it now appears that it's not so applicable to your problem, I think you can still refer to it. $\endgroup$
    – xzczd ♦
    Commented May 6, 2013 at 4:34
  • $\begingroup$ It's not surprising that the original initial condition becomes improper when the domain is changed, I guess the key point for solving your problem is to find a proper initial condition, which requires a lot of patience. $\endgroup$
    – xzczd ♦
    Commented May 6, 2013 at 7:05
  • $\begingroup$ @xzczd: Now I'm studying my problem better in order to overtake my difficulties. Thanks for all :) $\endgroup$
    – Petrus
    Commented May 7, 2013 at 1:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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