System of ODEs causes NDSolve::nlnum1 warning

Question

I am trying to solve a system of ODEs using some amplitude and frequency with NDSolve, but I get some warnings which is included here.

I think my program is fine until before SolAv. But, I could not find why this warning appears. The code I used is included here

Clear[LΛ, f, c, cn, A, v, An, As, vs, Q0]
(*  *)
LΛ[f_] := (8 π ν ρ f)/ρ;
c = cn An/A[x, t];
Eq1 = Simplify[D[A[x, t], t] + D[A[x, t] v[x, t], x]];

Eq2 = Simplify[
   A[x, t] D[v[x, t], t] + A[x, t] v[x, t] D[v[x, t], x] + 
     c^2 D[A[x, t], x] + LΛ[v[x, t]] /. 
    R -> Sqrt[A[x, t]/π]];


A[x_, t_] = As[x]; v[x_, t_] = vs[x];
vs[x_] = Q0/As[x];
LE = Simplify[{Eq1, Eq2}];
Sols = DSolve[{LE[[2]] == 0, As[L] == An}, As[x], x];
Simplify[Sols]
Am = {1.1290263207533879`*^-7 - 7.473554942563208`*^-8 I, 
   3.9210392313242154`*^-7 - 
    5.198540412010037`*^-7 I, -3.666731498809855`*^-8 + 
    4.026529625379335`*^-7 I, 
   2.907927596300974`*^-7 + 
    6.688553949022083`*^-7 I, -5.901727476937829`*^-7 - 
    6.780384037260203`*^-7 I, -3.4029233536625084`*^-7 - 
    2.1769713736464412`*^-7 I, 
   5.06797091099388`*^-7 - 
    2.5532173494757542`*^-8 I, -5.4970528551404145`*^-9 - 
    7.182784241041546`*^-8 I, -3.491221569640317`*^-7 + 
    2.6310572644437936`*^-7 I, -1.1440129576533528`*^-8 - 
    1.7690933852144399`*^-7 I};

ffp = {2/7, 4/7, 6/7, 8/7, 10/7, 12/7, 2, 16/7, 18/7, 20/7};
Clear[LΛ, f, c, cn, A, v, An, As, vs, Asf, vsf, Q0, \
ω, R, ν, ρ, EqIC, Qap, δ, vao, Aap]
AmpQ1 = {};
AmpQ2 = {};
TU = {};
TV = {};
TP = {};
TQ = {};
Print["Start"]
Table[
  cn = 111 ;(*m/s*)
  R0 = 4 10^-3; (*m*)
  Q0 = 2.03 /1000/60; (*m^3/s*)
  μ = 0.0147  ; (*pa-s*)
  ρ = 1167;(*kg/m^3*)
  ν = μ/ρ;
  L = 10;
  δ = 0.01;
  An = π R0^2;
  (*Tp=0.2;ω=2 π 1/Tp;*)
  ω = 2 π  ffp[[f]];
  (*Print["Disharge fluctuation=",δ Q0];*)
  
  M10C =
   1 - 2 BesselJ[
       1, α I^(3/2)]/(α I^(3/2)
         BesselJ[0, α I^(3/2)] );
  Par[α_] = I π μ α^2 (M10C^-1 - 1);
  G[x_] = 
   Par[α] /. α -> Sqrt[A[x, t] ω/(π*ν)];
  (*LΛ[f_]:=(8 π ν ρ f)/ρ;*)
  τ[x_, t_] = G[x]   v[x, t]/ρ;
  
  (*LΛ[f_]:=(8 π ν ρ f)/ρ;*)
  
  c = (cn An)/A[x, t];
  Eq1 = Simplify[D[A[x, t], t] + D[A[x, t] v[x, t], x]];
  
  Eq2 = Simplify[
    A[x, t] D[v[x, t], t] + A[x, t] v[x, t] D[v[x, t], x] + 
     c^2 D[A[x, t], x] + τ[x, t](*LΛ[v[x,
    t]]/.R\[Rule]Sqrt[A[x,t]/π]*)];
  
  A[x_, t_] = As[x] + ϵ Asf[x] E^(I ω t);
  v[x_, t_] = vs[x] + ϵ vsf[x] E^(I ω t);
  LEq = Series[{Eq1, Eq2}, {ϵ, 0, 1}] // Normal;
  EqIC = Coefficient[LEq, ϵ] /. t -> 0;
  (* steady flow solution *)
  As[x_] = 
   An E^(-((8 π Q0 (L - x) ν)/(-An^2 cn^2 + Q0^2)));(*As[x_]=
  An \[ExponentialE]^((7 π Q0 (L-x))/(12500 (An^2 cn^2-Q0^2)));*)
  
  vs[x_] = Q0/As[x];
  
  solAv = 
   First@NDSolve[
     Evaluate[{EqIC[[1]] == 0, EqIC[[2]] == 0, 
       Asf[L] == (- cn/I ω)  Derivative[1][Asf][L] (*Asf[L]*), 
       Asf[0] vs[0] + As[0] vsf[0] == Am [[f]] }], {Asf[
       x], vsf[x]}, {x, 0, L}];
  
  , {f, 1, Length[ffp]}];
Print["finish"];

The code causes the warning

NDSolve::nlnum1: The function value {(0. +0. I)-444/7 I π Asf'[10],-1.12903*10^-7+7.47355*10^-8 I} is not a list of numbers with dimensions {2} when the arguments are {0. +0. I,0. +0. I,0. +0. I,0. +0. I}.

...

If you help me to resolve the problem, it will be helpful for me.

Please help me.

Answer 1

First of all, I'd like to emphasize again, you should put some effort in creating a minimal working example. Anyway, the underlying issue of the messy sample turns out to be worth discussing, so let me give an answer.

The sample can be boiled down to the following single line:

NDSolve[{g'[x] == 0, g[0] == 0, f'[x] + f[x] == 0, f'[1] == 0}, {f, g}, {x, 0, 1}]

NDSolve::nlnum1: The function value {0.,f'[1]} is not a list of numbers with dimensions {2} when the arguments are {0.,0.,0.,0.}.

(* Input returned *)

Further check suggests this warning is caused by "Shooting" method:

NDSolve[{g'[x] == 0, g[0] == 0, f'[x] + f[x] == 0, f'[1] == 1234}, {f, g}, {x, 0, 1}, 
        Method -> {"Shooting", "StartingInitialConditions" -> {g[1] == 2345}}]

NDSolve::nlnum1: The function value {2345., - 1234 + f'[1]} is not a list of numbers with dimensions {2} when the arguments are {2.71828, 2345., 1., 2345.}.

(* Input returned *)

It's not clear to me how "Shooting" method parses the system, but it's clear to me "Shooting" method is having difficulty in handling the boundary condition f'[1] == 1234.

This isn't too surprising. Clearly, b.c. involving 1st order derivatvie is unusual for 1st order ODE: in cases like {f'[x] == x, f'[1] == 0} the solution won't even be uniquely determined!

"OK, but in our case the system happens to be well-determined, how to make NDSolve correctly parse the b.c.?" One possible way is to transform the unusual b.c. to a usual one: from f'[x] + f[x] == 0, f'[1] == 0 we can easily deduce f[1] == 0, so:

NDSolve[{g'[x] == 0, g[0] == 0, f'[x] + f[x] == 0, f[1] == 0}, {f, g}, {x, 0, 1}]

And this method is perfectly suitable for your original problem:

…
newbc = Eliminate[{EqIC[[1]] == 0, EqIC[[2]] == 0, 
                   Asf[L] == (-cn/I ω) Asf'[L]} /. x -> L, 
                  {Asf'[L], vsf'[L]}];

solAv = First@
  NDSolve[{EqIC[[1]] == 0, EqIC[[2]] == 0, newbc, 
     Asf[0] vs[0] + As[0] vsf[0] == Am[[f]]}, {Asf[x], vsf[x]}, {x, 0, L}]
…