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.
w^\[Prime]\[Prime])[z]
in "Out[32]" is wrong syntax. $\endgroup$Ra = 10^3; k = 1; \[Lambda] = 1;
. Second, replace(w^\[Prime]\[Prime])[z]
withw''[z]
. Third, add boundary condition forw[z]
, for example,w'[1/2]=0
. Finally, we have messages NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations. NDSolve::mconly: For the method NDSolve`IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions. $\endgroup$NDSolve
can't invert matrix and uses DAE instead. But DAE method is not applicable for complex system. Therefore, we can't solve this system withNDSolve
. $\endgroup$w^\[Prime]\[Prime])[z]
is caused by improper copy&paste, pasting it back to the notebook and press Ctrl+Shift+N, then the code will again become valid code. $\endgroup$