Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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;
\[Sigma]M = 0.00592251
\[Lambda]1 = 1.025;
\[Lambda]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] == \[Sigma]M,
y[0] == \[Lambda]1,
y[15] == \[Lambda]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] == \[Sigma]M,
      y[0] == \[Lambda]1,
      y[15] == \[Lambda]2}, y, {x, 0, 15}, 
      Method -> {"Shooting","StartingInitialConditions" -> {y[0] == \[Lambda]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?

share|improve this question

2 Answers 2

up vote 2 down vote accepted

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).

share|improve this answer
    
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? –  Petrus May 5 '13 at 14:09
    
@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. –  xzczd May 5 '13 at 16:21
    
@xzxzd thanks for your help –  Petrus May 5 '13 at 20:53
    
@Petrus See my edit. –  xzczd May 6 '13 at 4:49
    
xzczd:how the use of Table can help you? –  Petrus May 7 '13 at 1:42

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;
\[Sigma]M = 0.00592251;
\[Lambda]1 = 1.0253896074561006`;
\[Lambda]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] - \[Sigma]M;

LA = 0;
LB = 15;

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

Show[
     Plot[y[x] /. sol, {x, 0, LB}, PlotRange -> {Automatic, {1, 1.32}}, Frame -> True],
     Graphics[{
     Red,
     Point[{{LA, \[Lambda]1}, {LB, \[Lambda]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?

share|improve this answer
    
@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. –  xzczd May 5 '13 at 16:18
    
@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 –  Petrus May 5 '13 at 20:48
    
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. –  xzczd May 6 '13 at 4:34
    
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. –  xzczd May 6 '13 at 7:05
    
@xzczd: Now I'm studying my problem better in order to overtake my difficulties. Thanks for all :) –  Petrus May 7 '13 at 1:41

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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