Unable to find initial conditions that satisfy the residual function within specified tolerances

Question

I am trying to solve the differential system :

eqns = Flatten[{N[p5], N[p6], N[p7], N[p8]}]
(*
{k^2 T[z] + λ (-1. k^2 u[z] + du'[z]) - 
   0.842615 Sqrt[1/Ra] (k^4 u[z] + dddu'[z] - 2. k^2 du'[z]) + (
   17.4421 (-1. k^2 u[z] + k^2 w[z] + du'[z] - 1. w''[z]))/Sqrt[Ra] == 0., 
 ddu'[z] == dddu[z], du'[z] == ddu[z], 
 u'[z] == du[z], λ w[z] == (218.027 (u[z] - 1. w[z]))/Sqrt[Ra] + 
   84.2615 Sqrt[1/Ra] w'[z], λ T[z] + (18.0529 (T[z] - 1. Tp[z]))/Sqrt[
   Ra] + (-(0.0863887/2.71828^(4.22949 (0.5 - 1. z))) - 
      0.101592 2.71828^(3.59656 (0.5 - 1. z))) u[z] - (
   1.18678 (-1. k^2 T[z] + dT'[z]))/Sqrt[Ra] == 0., 
 T'[z] == dT[
   z], -((53.3319 (T[z] - 1. Tp[z]))/Sqrt[
    Ra]) + λ Tp[
     z] + (0.015203/2.71828^(4.22949 (0.5 - 1. z)) - 
      0.015203 2.71828^(3.59656 (0.5 - 1. z))) w[z] - 84.2615 Sqrt[1/Ra] Tp'[z] == 
  0.}
*)

with boundary conditions :

ics = Join[p1, p3] /. {cte1 -> 1 + I*1, cte2 -> 1 + I*1}
(* {u[1/2] == 0, du[1/2] == 0, ddu[1/2] == 1 + I, 
 dddu[1/2] == 1 + I, T[1/2] == 0, dT[1/2] == 1, w[1/2] == 0, 
 Tp[1/2] == 0} *)

And I get the following error message:

NDSolve[
 Join[eqns, ics], {u, du, ddu, dddu, w, T, dT, Tp}, {z, -1/2, 1/2}]

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.

NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.

(* {} *)

Can someone please help me solving this problem?

Thank you.

Answer 1

As pointed out by Alex, the parameters are missing. Alex kindly tried guessing the missing values, but parameters sometimes can severely influence the solving of equation system, so it's better to show us your specific parameters. Anyway, aside from the parameter issue, the system shows something interesting, so let me post an answer using parameters in Alex's answer.

As already shown in question, NDSolve is having difficulty in automatically calling an ODE solver, which seems to be strange (Well, perhaps it's not that strange, we already have example like this), so I decide to try manually eliminating the variables in the system. Luckily the elimination is simple enough for this system:

Clear[k, Ra, λ]
dddu[z_] = ddu'[z];
ddu[z_] = du'[z];
du[z_] = u'[z];
dT[z_] = T'[z];

neweq = {k^2 T[z] + λ (-1. k^2 u[z] + du'[z]) - 
     0.842615 Sqrt[1/Ra] (k^4 u[z] + dddu'[z] - 2. k^2 du'[z]) + (
     17.4421 (-1. k^2 u[z] + k^2 w[z] + du'[z] - 1. w''[z]))/Sqrt[Ra] == 0.,
   λ w[z] == (218.027 (u[z] - 1. w[z]))/Sqrt[Ra] + 84.2615 Sqrt[1/Ra] w'[z],
   λ T[z] + (18.0529 (T[z] - 1. Tp[z]))/Sqrt[
     Ra] + (-(0.0863887/2.71828^(4.22949 (0.5 - 1. z))) - 
        0.101592 2.71828^(3.59656 (0.5 - 1. z))) u[z] - (
     1.18678 (-1. k^2 T[z] + dT'[z]))/Sqrt[Ra] == 0.,
   -((53.3319 (T[z] - 1. Tp[z]))/Sqrt[
      Ra]) + λ Tp[
       z] + (0.015203/2.71828^(4.22949 (0.5 - 1. z)) - 
        0.015203 2.71828^(3.59656 (0.5 - 1. z))) w[z] - 
        84.2615 Sqrt[1/Ra] Tp'[z] == 0.};

ic = {u[1/2] == 0, du[1/2] == 0, ddu[1/2] == 1 + I, dddu[1/2] == 1 + I, 
      T[1/2] == 0, dT[1/2] == 1, w[1/2] == 0(*,w'[1/2]==0*), Tp[1/2] == 0};

rule = Solve[
     Most@neweq, {u, T, Tp}[z] // Through][[1]] /. (h_[z] -> rhs_) :> (h -> 
      Function @@ {z, rhs});

neweqfinal = Last@neweq //. rule;

icfinal = ic //. rule // Simplify[#, {Ra > 0, k > 0}] &;

Let's have a quick check:

Cases[neweqfinal, _[z], Infinity] // Union
(* {w[z], w'[z], w''[z], w'''[z], w''''[z], Derivative[5][w][z], 
 Derivative[6][w][z], Derivative[7][w][z], Derivative[8][w][z]} *)

Cases[icfinal, _[1/2], Infinity] // Union
(* {w[1/2], w'[1/2], w''[1/2], w'''[1/2], w''''[1/2], Derivative[5][w][1/2], 
 Derivative[6][w][1/2], Derivative[7][w][1/2]} *)

So it seems that, after elimination we obtain a eighth order ODE with 8 i.c.s, NDSolve should be able to handle it.

And it can indeed handle it:

Ra = 10^3; k = 1; λ = 1;

solw = NDSolveValue[{neweqfinal, icfinal}, w, {z, -1/2, 1/2}, 
    InterpolationOrder -> All]; // AbsoluteTiming

Clear[solu, solT, solTp]

{solu[z_], solT[z_], solTp[z_]} = {u[z], T[z], Tp[z]} //. rule /. w -> solw;

ReImPlot[{solu[z], solT[z], solTp[z], solw[z]}, {z, -1/2, 1/2}]

The solution looks different from Alex's, I'm not immediately sure about reason.

---

Starting from the system used by Alex, we can obtain the same result using the method above. Notice w'[1/2]==0 is removed because it's redundant in this method:

Clear[du, ddu, dddu, dT, Ra, k, λ]
eqs = Simplify[{w'[z] == w1[z], 
     k^2 T[z] + λ (-1. k^2 u[z] + du'[z]) - 
       0.842615 Sqrt[1/Ra] (k^4 u[z] + dddu'[z] - 2. k^2 du'[z]) + (
       17.4421 (-1. k^2 u[z] + k^2 w[z] + du'[z] - 1. w1'[z]))/Sqrt[Ra] == 0., 
     ddu'[z] == dddu[z], du'[z] == ddu[z], 
     u'[z] == du[z], λ w[z] == (218.027 (u[z] - 1. w[z]))/Sqrt[Ra] + 
       84.2615 Sqrt[1/Ra] w'[z], λ T[z] + (18.0529 (T[z] - 1. Tp[z]))/Sqrt[
       Ra] + (-0.0863887 2.71828^(-4.22949 (0.5 - 1. z)) - 
          0.101592 2.71828^(3.59656 (0.5 - 1. z))) u[z] - (
       1.18678 (-1. k^2 T[z] + dT'[z]))/Sqrt[Ra] == 0., 
     T'[z] == dT[
       z], -((53.3319 (T[z] - 1. Tp[z]))/Sqrt[
        Ra]) + λ Tp[
         z] + (0.015203 2.71828^(-4.22949 (0.5 - 1. z)) - 
          0.015203 2.71828^(3.59656 (0.5 - 1. z))) w[z] - 
       84.2615 Sqrt[1/Ra] Tp'[z] == 0.}] // Rationalize[#1, 0] &;
ics = {u[1/2] == 0, du[1/2] == 0, ddu[1/2] == 1 + I, dddu[1/2] == 1 + I, 
       T[1/2] == 0,dT[1/2] == 1, w[1/2] == 0, Tp[1/2] == 0};    
rule2 = Solve[
     Most@eqs, {u, T, Tp, du, ddu, dddu, dT, w1}[z] // Through][[1]] /. (h_[z] -> 
      rhs_) :> (h -> Function @@ {z, rhs});   
neweqfinal2 = Last@eqs //. rule2;   
icfinal2 = ics //. rule2 // Simplify[#, {Ra > 0, k > 0}] &;
Ra = 10^3; k = 1; λ = 1;   
solw2 = NDSolveValue[{neweqfinal2, icfinal2}, w, {z, -1/2, 1/2}, 
    InterpolationOrder -> All, WorkingPrecision -> 16]; // AbsoluteTiming    
{solu2[z_], solT2[z_], solTp2[z_]} = {u[z], T[z], Tp[z]} //. rule2 /. w -> solw2;    
ReImPlot[{solu2[z], solT2[z], solTp2[z], solw2[z]}, {z, -1/2, 1/2}]
ics //. rule2 /. w -> solw2
(* {True, True, True, True, True, True, True, True} *)

---

Aha, after thinking carefully about why NDSolve is having difficulty in automatically calling an ODE solver, I found another approach. Let's take a closer look at the system eqs. Clearly it's linear system involving 9 different 1st order derivative term, but:

Union@Cases[#, HoldPattern@Derivative[__][__][__], Infinity] & /@ eqs
(* {{w'[z]}, 
    {dddu'[z], du'[z], w1'[z]}, 
    {ddu'[z]}, 
    {du'[z]}, 
    {u'[z]}, 
    {w'[z]}, 
    {dT'[z]}, 
    {T'[z]}, 
    {Tp'[z]}} *)

1st and 6th equation only involves w'[z] term, w1'[z] term and dddu'[z] term only exist in 2nd equation, so it'll be impossible to transform the system to the "standard form" required by an ODE IVP solver only by algebraic elimination. (See this post to learn more about the "standard form". ) This is probably the reason why NDSolve is having difficulty in automatic transformation.

Then can we do something to help NDSolve? The answer is yes, we can easily generate another equation involving derivative of w1 from 1st and 6th equation in the following manner:

eqmid = Eliminate[eqs[[{1, 6}]], w'[z]];
eqadd = D[eqmid, z]
icadd = eqmid /. z -> 1/2 /. Rule @@@ ics // Simplify
(* 168523 Sqrt[10]
   w1'[z] == -2 (218027 Sqrt[10] u'[z] - 100000 w'[z] - 218027 Sqrt[10] w'[z]) *)
(* w1[1/2] == 0 *)

Now we can use NDSolve to solve the problem as follows:

{solu3, solT3, solTp3, solw3} = 
   NDSolveValue[{Rest@eqs, eqadd, ics, icadd}, {u, T, Tp, w}, {z, -1/2, 
     1/2}]; // AbsoluteTiming

ReImPlot[{solu3[z], solT3[z], solTp3[z], solw3[z]}, {z, -1/2, 1/2}]

Definitions of eqs and ics are the same as above. The obtained ReImPlot is the same as above.

Answer 2

The matrix of this system is singular, it is why it can't be solved with NDSolve directly since DAE method is applicable to the machine real code only. To solve this system, we can use Euler wavelets collocation method as follows. First, we turn this system to the first order by definition a new function w1[z]==w'[z] with initial condition w1[.5]==0. Then we have (this is numerical example only)

eqs = {w'[z] == w1[z], 
    k^2 T[z] + λ (-1. k^2 u[z] + Derivative[1][du][z]) - 
      0.842615 Sqrt[
        1/Ra] (k^4 u[z] + Derivative[1][dddu][z] - 
         2. k^2 Derivative[1][du][z]) + (17.4421 (-1. k^2 u[z] + 
           k^2 w[z] + Derivative[1][du][z] - 1. w1'[z]))/Sqrt[Ra] == 
     0., Derivative[1][ddu][z] == dddu[z], 
    Derivative[1][du][z] == ddu[z], 
    Derivative[1][u][z] == 
     du[z], λ w[z] == (218.027 (u[z] - 1. w[z]))/Sqrt[Ra] + 
      84.2615 Sqrt[1/Ra] Derivative[1][w][z], λ T[
        z] + (18.0529 (T[z] - 1. Tp[z]))/
       Sqrt[Ra] + (-0.0863887 2.71828^(-4.22949 (0.5 - 1. z)) - 
         0.101592 2.71828^(3.59656 (0.5 - 1. z))) u[
        z] - (1.18678 (-1. k^2 T[z] + Derivative[1][dT][z]))/
       Sqrt[Ra] == 0., 
    Derivative[1][T][z] == 
     dT[z], -((53.3319 (T[z] - 1. Tp[z]))/Sqrt[Ra]) + λ Tp[
        z] + (0.015203 2.71828^(-4.22949 (0.5 - 1. z)) - 
         0.015203 2.71828^(3.59656 (0.5 - 1. z))) w[z] - 
      84.2615 Sqrt[1/Ra] Derivative[1][Tp][z] == 0.} // Simplify;
bcs = {u[1/2] == 0, du[1/2] == 0, ddu[1/2] == 1 + I, 

dddu[1/2] == 1 + I, w[1/2] == 0, w1[1/2] == 0, T[1/2] == 0, 
   dT[1/2] == 1, Tp[1/2] == 0};

Ra = 10^3; k = 1; λ = 1; 

Second, we compute matrix

var = {u'[z], du'[z], ddu'[z], dddu'[z], w'[z], w1'[z], T'[z], dT'[z],
   Tp'[z]}
var0 = {u[z], du[z], ddu[z], dddu[z], w[z], w1[z], T[z], dT[z], Tp[z]};
UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] := 
 Piecewise[{{2^(k/2) UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t < 
     n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 4; With[{k = k0, M = M0}, 
 nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}]; xcol = 
 Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
M = nn; A = Array[a, {9, M}]; B = Array[b, {9}];
eqn = Flatten[
  Table[eqs /. 
     Flatten[Table[{var[[i]] -> A[[i]] . Psi[xcol[[j]]], 
        var0[[i]] -> A[[i]] . int1[xcol[[j]]] + B[[i]]}, {i, 
        9}]] /. {z -> xcol[[j]] - 1/2}, {j, M}]]; bc = 
 Table[A[[i]] . int1[1.] + B[[i]] == bcs[[i, 2]], {i, 9}];
varM = Join[Flatten[A], B]; {vec, mat} = 
 CoefficientArrays[Join[eqn, bc], varM];

Solution

sol = LinearSolve[mat, -vec];rul = Table[varM[[i]] -> sol[[i]], {i, Length[varM]}];

Visualization

 ReImPlot[
  Evaluate[
   Table[(A[[i]] . int1[z + .5] + B[[i]]) /. 
     rul, {i, {1, 5, 7, 9}}]], {z, -1/2, 1/2}]

Note, this solution exactly the same as computed by xzczd. Thanks to him we know solution, and my first attempt with PseudoInverse has an error.