1
$\begingroup$

I would greatly appreciate your help in implementation of Mathematica for solving a 2nd order differential equation expressing some important engineering problem.

Apparently, this equation has no analytical solution. Therefore I attempt to solve it numerically with maximal possible approach to the analytical effect, while using the ParametricNDSolve command.

The equation has 3 variables: n,Ω and ψ and differentiation is by z ($0<z<1$). n can be above zero and below it, however the values rarely exceed ±3. Ω and ψ are positive numbers and they can reach big values. As the 1st approach I have defined n=3. The expression is:

Solution=ParametricNDSolve[{Ω*f''[z][z]-f'[z]-ψ*f[z]^n==0,f[0]-Ω*f'[0]==1,f'[1]==0},f,{z,0,1},{Ω,ψ}];

Unfortunately Mathematica is unable to solve this problem in wide ranges of Ω and ψ and produces the following messages:

ParametricNDSolve::ndsz: At z$349 == 0.9049990200792034`, step size is effectively zero; singularity or stiff system suspected. >>
ParametricNDSolve::ndsz: At z$343 == 0.8954517351902076`, step size is effectively zero; singularity or stiff system suspected. >>
ParametricNDSolve::ndsz: At z$343 == 0.8548633486473108`, step size is effectively zero; singularity or stiff system suspected. >>
General::stop: Further output of ParametricNDSolve::ndsz will be suppressed during this calculation. >>

And the results at some points are non-physical.

Actually, the practically important solution is at z=1. I have tried that in the command: {z,1,1} and did not get these messages. However they have appeared while attempting to create the 3 dimensional plot:

Plot3D[Evaluate[f[Ω,ψ][1]/.Solution],{Ω,1/5,20},{ψ,0,7}]

And no progress has been achieved.

After reading answers to other questions on the command ParametricNDSolve I have tried to define a solution method:

Method→"StiffnessSwitching","ExtrapolationHandler"→{Indeterminate&,"WarningMessage"→False}

Unfortunately it did not result in any progress as well.

I would greatly appreciate your opinions if something else can be done in this case in order to overcome the singularity. Thank you a lot in advance.

$\endgroup$
  • $\begingroup$ There are syntax errors in your code, such as \[Psi]*f[z] instead of ψf[z], \[CapitalOmega]*f''[z] instead of Ωf"[z] and f''[z] instead of f"[z]. The latter means that '' this is twice Prim, but not the quote sign. After I corrected all this, your equation nicely solves. I propose to close this question as the off-topic. The problem is due to basic syntax errors. $\endgroup$ – Alexei Boulbitch Sep 28 '15 at 9:50
  • $\begingroup$ Dear Alexei, thank you a lot for your comment. Unfortunately, the expression just has not been correctly reproduced by copying from my Mathematica code and I have edited my question. In fact there are no syntax errors there, and the procedure runs, yielding correct results at some points. Unfortunately, there is a problem of singularity at other points, where the results are non-physical. $\endgroup$ – Sitra Ahra Sep 28 '15 at 9:59
  • $\begingroup$ @SitraAhra For copying & formatting code, you may find this this meta Q&A helpful. $\endgroup$ – Michael E2 Sep 28 '15 at 10:32
  • $\begingroup$ Dear Michael, thank you a lot for this useful comment. I will implement your important suggestion in my next questions. $\endgroup$ – Sitra Ahra Sep 28 '15 at 10:44
1
$\begingroup$

Fast and dirty, try this:

 Clear[\[CapitalOmega], \[Psi], sol];
n = 3;
\[Psi] = 1;
lst = Flatten[
  Table[sol = 
     NDSolve[{\[CapitalOmega]*f''[z] - f'[z] - \[Psi]*f[z]^n == 0, 
       f[0] - \[CapitalOmega]*f'[0] == 1, f'[1] == 0}, f, {z, 0, 1}];
    Table[
     {\[CapitalOmega], z, f[z] /. sol}, {z, 0, 1, 
      0.1}], {\[CapitalOmega], 1, 2, 0.1}] /. {x_, y_, {z_}} -> {x, y,
      z}, 1]

Now one may plot the obtained list, lst:

    ListPlot3D[lst, 
 AxesLabel -> {Style["\[CapitalOmega]", 18, Italic], 
   Style["z", 18, Italic], Style["Res", 18, Italic]}]

yielding this:

enter image description here

Have fun!

$\endgroup$
  • $\begingroup$ Dear Alexei, thank you a lot! Unfortunately this procedure does not solve the singularity problem. For example taking [Psi] = 4 (and bigger) the same error messages appear and some results are again non-physical. At the same time the problem could be simplified by considering z=1 only. The problem is that perhaps I do not know how to do it. Or perhaps they are some addtional methods. Thank you a lot again. $\endgroup$ – Sitra Ahra Sep 28 '15 at 13:39
  • $\begingroup$ @Sitra Ahra Well, it seems that this equation, indeed, becomes singular at \[Psi]=3.92... Nothing to do. I tried the Method -> "StiffnessSwitching" it does not help. $\endgroup$ – Alexei Boulbitch Sep 28 '15 at 14:36

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.