3
$\begingroup$

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]
$\endgroup$
  • $\begingroup$ It is strange that NDSolve stumbles on a simple system of equations. $\endgroup$ – Alex Trounev Oct 8 '19 at 10:28
  • $\begingroup$ 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. $\endgroup$ – sixpenny Oct 8 '19 at 10:59
  • $\begingroup$ How about my first question to solve T[t,z,r] and S[t,z,r]? Thank you. $\endgroup$ – sixpenny Oct 8 '19 at 11:29
  • $\begingroup$ NDSolve[{pde, ic}, {S, T}, {t, 0, 10}, {z, 0, 10}, DependentVariables -> {T[t, z], S[t, z]}] $\endgroup$ – user21 Oct 8 '19 at 11:47
  • $\begingroup$ With the new added equation, the first system is over-determined. $\endgroup$ – xzczd Oct 8 '19 at 11:48
3
$\begingroup$

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,
   please find it in the link above. *)
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}]

enter image description here

$\endgroup$
1
$\begingroup$

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.

addendum(might be the answer)

You can solve the problem numerically in two steps:

solS = NDSolveValue[{D[S[t, z, r], z] == -S[t, z, r],S[t, 0, r] == Exp[-r^2]}, S, {t, 0, 10}, {z, 0, 10}, {r, 0, 10} ,Method -> "FiniteElement" ]
solT = NDSolveValue[{D[T[t, z, r], t] == 1/n solS[t, z, r], T[0, z, r] == 1}, T, {t, 0, 10}, {z, 0, 10}, {r, 0, 10} ] 
ContourPlot[solT[t, z, 1], {z, 0, 10}, {t, 0, 10},PlotLegends -> Automatic]

enter image description here

$\endgroup$
  • $\begingroup$ Since the system is linear, this should be the only analytic solution. I've found it with LaplaceTransform and DSolve, see my answer for more detail. $\endgroup$ – xzczd Oct 8 '19 at 14:17
  • $\begingroup$ @xzczd, I will read your great code carefully, thanks a lot. Never thought that I should write so many codes in mma. $\endgroup$ – sixpenny Oct 8 '19 at 14:20
  • $\begingroup$ @Ulrich Neumann, I will also try your way. In fact, I simplied my pdes. The real equation is D[S[t, z, r], z] == - 2 Im[(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)]* S[t, z, r] and D[T[t, z, r], t] == -n^-1*D[S[t, z, r], z]. I wonder how to solve it in your way. $\endgroup$ – sixpenny Oct 8 '19 at 15:18
  • $\begingroup$ @sixpenny This system is completely different to that you asked for! $\endgroup$ – Ulrich Neumann Oct 8 '19 at 15:23
  • $\begingroup$ @UlrichNeumann, Sorry. At the beginning, I wanted to simply the pdes only to know the skills to solve it. Did not know that they are so different. The real pdes is updated in the original post. Anyway, I also enjoy your skill for the simplified pdes. Thanks. $\endgroup$ – sixpenny Oct 8 '19 at 15:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.