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. Actually, 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] - f'[z] - ψ f[z]^n == 0,
                              f[0] - Ω f'[0] == 1,f'[1] == 0}, f, {z, 1, 1}, {Ω, ψ}];

Unfortunately Mathematica is unable restrict itself to z = 1 and produces the singularity messages at z < 1.

My question is thus following: is it possible to restrict the solution to z=1 without attempting to solve at z<1?

  • $\begingroup$ Is it f''[z][z] or f''[z] or z f''[z] ?? $\endgroup$
    – thils
    Nov 9 '15 at 10:44
  • $\begingroup$ it is f''[z], sorry for this typo $\endgroup$
    – Sitra Ahra
    Nov 9 '15 at 11:05
  • $\begingroup$ Because one boundary condition is at z = 0 and the other at z = 1, ParametricNDSolve must integrate from 0 to 1. $\endgroup$
    – bbgodfrey
    Nov 10 '15 at 3:26

Because boundary conditions are specified at z = 0 and at z = 1, ParametricNDSolve must integrate between those two endpoints. However, this does not in general prevent obtaining a solution.

s = ParametricNDSolveValue[{Ω f''[z] - f'[z] - ψ f[z]^3 == 0, 
    f[0] - Ω f'[0] == 1, f'[1] == 0}, f, {z, 0, 1}, {Ω, ψ}];

which then can be evaluated at s[Ω, ψ][1]. For instance,

Plot[Evaluate[Table[s[Ω, ψ][1], {ψ, 5}]], {Ω, 1, 5}, AxesLabel -> {Ω, f}]

enter image description here

It is the case that stiffness is encountered at small Ω for ψ == 5.


The parameter ψ can be eliminated by rescaling f by ψ^(-1/(n - 1)), which reduces parameter space from two to one dimensions. With the rescaled f designated g,

t = ParametricNDSolveValue[{Ω g''[z] - g'[z] - g[z]^3 == 0, g[0] - Ω g'[0] == 1, 
    g'[1] == 0}, g, {z, 0, 1}, {Ω}]
Plot[t[Ω][1], {Ω, 1, 20}, PlotRange -> {{0, 20}, All}, AxesLabel -> 
    {"Ω", "f \!\(\*SuperscriptBox[\(ψ\), FractionBox[\(1\), \(n - 1\)]]\)"

enter image description here

  • $\begingroup$ Thank you a lot! I have the following question: perhaps the singularity could be avoided or reduced by defining other solution method? $\endgroup$
    – Sitra Ahra
    Nov 10 '15 at 8:03
  • $\begingroup$ @SitraAhra Although there probably are several ways to avoid this issue, I recommend the rescaling shown in the addendum above. $\endgroup$
    – bbgodfrey
    Nov 11 '15 at 2:53
  • $\begingroup$ Thank you a lot! Indeed the rescaling solves the problem, for all positive values of Ω and n. Thank you a lot again! $\endgroup$
    – Sitra Ahra
    Nov 13 '15 at 8:52

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.