# Need more equations in PDEs by NDSolve? Bug of mma?

I am trying to solve set of pdes as below

ClearAll["Global*"] ;
n = 100;
pde = {D[S[t, z, r], z] == -S[t, z, r],
D[T[t, z, r], t] == n^-1 S[t, z, r]};
ic = {T[0, z, r] == 1, S[t, 0, r] == Exp[-r^2]};
sol = NDSolve[{pde, ic}, {S, T}, {t, 0, 10}, {z, 0, 10}, {r, 0, 10}];
ContourPlot[Evaluate[T[t, z, 1] /. sol], {z, 0, 10}, {t, 0, 10}, PlotLegends -> Automatic]

Do I need more equations and initial conditions? By the way, the pdes can not be solved when sequence of equation is different. Is it bug of mma? For example, pdes below can not be solved.

ClearAll["Global*"] ;
n = 100;
pde = {D[T[t, z], t] == n^-1 S[t, z], D[S[t, z], z] == -S[t, z]};
ic = {T[0, z] == 1, S[t, 0] == 1};
sol = NDSolve[{pde, ic}, {S, T}, {t, 0, 10}, {z, 0, 10}];
ContourPlot[Evaluate[T[t, z] /. sol], {z, 0, 10}, {t, 0, 10}, PlotLegends -> Automatic]

while pdes below can be solved.

ClearAll["Global*"] ;
n = 100;
pde = {D[S[t, z], z] == -S[t, z], D[T[t, z], t] == n^-1 S[t, z]};
ic = {T[0, z] == 1, S[t, 0] == 1};
sol = NDSolve[{pde, ic}, {S, T}, {t, 0, 10}, {z, 0, 10}];
ContourPlot[Evaluate[T[t, z] /. sol], {z, 0, 10}, {t, 0, 10}, PlotLegends -> Automatic]

In fact, my real pdes is a little more complicated than 1st pde.

ClearAll["Global*"] ;
n = 100;B=10;p=10^-8;T0=10^-5;S0=10^-6;
pde = {D[S[t, z, r], z] == -2Im[k]*S[t, z, r],
D[T[t, z, r], t] == n^-1 S[t, z, r]}/.k->(1 - (n*T[t, z, r]^(3/2))/(T[t, z, r]^(3/2) - B*T[t, z, r]^(3/2) + I*p*n))^(1/2);
ic = {T[0, z, r] == T0, S[t, 0, r] == S0*Exp[-r^2]};
sol = NDSolve[{pde, ic}, {S, T}, {t, 0, 10}, {z, 0, 10}, {r, 0, 10}];
ContourPlot[Evaluate[T[t, z, 1] /. sol], {z, 0, 10}, {t, 0, 10}, PlotLegends -> Automatic]
• It is strange that NDSolve stumbles on a simple system of equations. Oct 8, 2019 at 10:28
• I just asked experts in PDE. They said that it is better to write code by myself to solve pdes, not using mma or library in matlab. In this way, I know exactly what is happening. Oct 8, 2019 at 10:59
• How about my first question to solve T[t,z,r] and S[t,z,r]? Thank you. Oct 8, 2019 at 11:29
• NDSolve[{pde, ic}, {S, T}, {t, 0, 10}, {z, 0, 10}, DependentVariables -> {T[t, z], S[t, z]}] Oct 8, 2019 at 11:47
• With the new added equation, the first system is over-determined. Oct 8, 2019 at 11:48

Since your target is just to solve 1st PDE system numerically, let me show you a finite difference method (FDM) based solution that should be applicable for your earlier problems (1, 2), as long as they're correct and well-posed.

First we solve your toy example analytically for comparision:

n = 100;
pde = {D[S[t, z, r], z] == -S[t, z, r], D[T[t, z, r], t] == n^-1 S[t, z, r]};
ic = {T[0, z, r] == 1, S[t, 0, r] == Exp[-r^2]};

rule = HoldPattern@LaplaceTransform[a_, __] :> a;
{teq, tic} = LaplaceTransform[{pde, ic[[2]]}, t, s] /. Rule @@ ic[[1]] /. rule;
{asolS[t_, z_, r_], asolT[t_, z_, r_]} =
InverseLaplaceTransform[DSolve[{teq, tic}, {S[t, z, r], T[t, z, r]}, {z}][[1, All, -1]],
s, t]
(* {E^(-r^2 - z), 1/100 E^(-r^2 - z) (100 E^(r^2 + z) + t)} *)
{pde, ic} /. {S -> asolS, T -> asolT}
(* {{True, True}, {True, True}} *)

Then we solve the problem based on FDM. I'll use pdetoae for the generation of finite difference equation:

points = 25; domain = {0, 10}; grid = Array[# &, points, domain]; difforder = 2;
(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[{S, T}[t, z, r], {grid, grid, grid}, difforder];
del = Rest;
ae = {del /@ ptoafunc@pde[[1]], del@ptoafunc@pde[[2]]};
aeic = ptoafunc@ic;
var = Outer[#[#2, #3, #4] &, {S, T}, grid, grid, grid, 1] // Flatten;

{barray, marray} = CoefficientArrays[{ae, aeic} // Flatten, var]

sollst = LinearSolve[marray, -N@barray];

{solS, solT} =
ListInterpolation[#, {grid, grid, grid}] & /@
ArrayReshape[sollst, {2, points, points, points}]

Manipulate[Plot[{#[t, z, r], #2[t, z, r]}, {z, 0, 10}, PlotRange -> All,
PlotStyle -> {Automatic, {Red, Dashed}}] & @@@ {{asolS, solS}, {asolT, solT}}, {t, 0,
1}, {r, 0, 1}]

An extended comment (might be a special analytical solution):

Your pde in S D[S[t, z, r], z] == -S[t, z, r] together with condition S[t, 0, r] == Exp[-r^2] has the general solution S[t,z,r]=Exp[-z] Exp[-r^2]

{D[S[t, z, r], z] == -S[t, z, r], S[t, 0, r] == Exp[-r^2]} /. S -> Function[{t, z, r}, Exp[-z] Exp[-r^2] ] // Simplify
(*{True,True}*)

That means S is independent of t!

The remaining equations give T[t,z,r]==1+t 1/n Exp[-z] Exp[-r^2]

{D[T[t, z, r], t] == n^-1 Exp[-z] Exp[-r^2], T[0, z, r] == 1} /. T -> Function[{t, z, r}, 1 + t/n Exp[-z] Exp[-r^2]] // Simplify
(*{True,True}*)

Don't know why Mathematica couldn't find this solution.