NDSolve with two parameters

Question

I was trying to solve a ODE numerically. It has two parameters (w and z0) which I want to vary. The following code gives an error: "Cannot find starting value for the variable".

g[z_] := 1 + 3 ( (z/z0))^4 - 4 ( (z/z0))^3 
DE = x''[z] + D[g[z]/z^2, z]/(g[z]/z^2) x'[z] +  (w^2)/(g[z])^2 x[z] == 0;

a = 10^-4;
zb = 0.0001;

pfun = ParametricNDSolveValue[{
            DE, x'[z0 - a] == - I (z0 w)/(6 w^2) x[z0 - a], x[zb] == 1
           }, x'[zb], {z, zb, z0 - a}, {w, z0}]
Plot[Evaluate[Table[
          Re[pfun[w, z0]]
    , {z0, 0.5, 1, .1}]]
    , {w,-0.01,0.01}, PlotRange -> All]

This error is an "Initial Condition" definition error. Please consider rewriting x'[z0 - a] == - I (z0 w)/(6 w^2) x[z0 - a], x[zb] == 1 in a more clear way. — Serhan Aya, Nov 10 '15 at 07:16
@SerhanAya Thanks for your comment! But I have no idea how to modify it. Can you suggest something which might work? — Physics Moron, Nov 10 '15 at 11:09

score 3 · Accepted Answer · answered Nov 10 '15 at 11:08

3

This seems to be caused by the boundary condition evaluated at a variable point. As an illustration, this

y[4][3] /. ParametricNDSolve[{
              y''[x] + y'[x] + y[x] == 0, 
              y'[10] == y[10], 
              y[0] == u
           }, y, {x, 0, 10}, u]

will evaluate just fine, but changing the y' initial condition to

y[4][3] /. ParametricNDSolve[{
              y''[x] + y'[x] + y[x] == 0, 
              y'[u] == y[u], 
              y[0] == u
           }, y, {x, 0, 10}, u]

will return the error

ParametricNDSolve::ndsv: Cannot find starting value for the variable y'.

In your case, of course, you've got your boundary conditions and that's that, and you can't really change them without changing your problem. What you can do is either

move away from ParametricNDSolve as a solver, in favour of functionalized constructs such as

f[u_?NumericQ] := f[u] = y /. First[NDSolve[{
                                       y''[x] + y'[x] + y[x] == 0,
                                       y'[u] == y[u], 
                                       y[0] == u
                                      }, y, {x, 0, 10}]]

which evaluates f[4][3] just fine, or

re-scale your problem to a new independent variable ζ with a fixed range, such as ζ=z/(z0-a-zb), for which this shouldn't be a problem.

answered Nov 10 '15 at 11:08

Emilio Pisanty

10,255
1
36
69

Thanks a lot for your answer! I have some very basic doubts. Why the "First" is so crucial here? Also if I have two such parameter how to plot the function as I did in my question? – Physics Moron Nov 10 '15 at 12:26
1

The First is because NDSolve returns its results as a list of lists of replacement rules (as a parallel to DSolve, which can return multiple independent solutions, cf. the output of DSolve[y'[x] + x y'[x]^2 == 1, y, x]), as opposed to ParametricNDSolve, which returns a simple list of replacement rules. If you have more than one parameter simply use f[u_?NumericQ,v_?NumericQ] := f[u,v] = (whatever). For more on _?NumericQ, see here. – Emilio Pisanty Nov 10 '15 at 13:21
Got it now. Another quick confusion : Plot[Table[f[u][x],{u,1,10}],{x,0,10}] gives back all curves but with same color. Any way out how to get different colors for different curves? – Physics Moron Nov 10 '15 at 13:39
1

Wrap the Table inside an Evaluate. This will force the Table to evaluate to an explicit list, which the Plot command then recognizes as things that need to be plotted in different styles. It's admittedly a bit of a bizarre idiom. – Emilio Pisanty Nov 10 '15 at 14:52
In my original question if I use the following code (following your suggestion), it doesn't give me the output. I must be making some silly mistake here!
f[w_?NumericQ, z0_?NumericQ] := f[w, z0] = x /. First[NDSolve[{DE, x'[z0-a] == -I (z0 w)/(6 a^2) x[z0-a], x[zb] == 1},x ,{z,zb,z0-a}]]

f[0.001, 2][zb]
– Physics Moron Nov 10 '15 at 19:29
Always test with a specific value before you run with it. Try Block[{w = 0.001, z0 = 2}, NDSolve[{DE, x'[z0 - a] == -I (z0 w)/(6 a^2) x[z0 - a], x[zb] == 1}, x, {z, zb, z0 - a}] ]. It returns the error NDSolve::bvluc: The equations derived from the boundary conditions are numerically ill-conditioned. The boundary conditions may not be sufficient to uniquely define a solution. This is a fundamental problem with your equation and you need to address it directly; there's little that we can do to help without knowing what it is you're working on. – Emilio Pisanty Nov 10 '15 at 19:48

NDSolve with two parameters

1 Answers1