The user-defined function pdetoode and pdetoae, developed by @xzczd, is very useful to deal with a PDE system when there is difficult to discretize it into a system of ODEs. In this situation, NDSolve often transforms it into a DAE system and then uses a DAE solver, which is weaker than ODE solver. That is why it normally failed in such a case.
Here I tried to use pdetoae to solve such a PDE system with random initial conditions. Then I want to use the solutions to solve for another equation further by calling the previous solutions somehow. Please find the function pdetoode and pdetoae at this link.
Clear[f, m, Tend];
f[y_] := y; m = 300; Tend = 10;
(*Ramdom Iintial Conditions*)
Clear[a0, b0, c0, d0]
SeedRandom[1];
a0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
SeedRandom[2];
b0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
SeedRandom[3];
c0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
SeedRandom[4];
d0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
Clear[difforder, domain, points, grid];
difforder = 4;
domain[y] = {-1, 1}; domain[t] = {0, Tend};
points[y] = 201; points[t] = 101;
(grid@# = Array[# &, points@#, domain@#]) & /@ {y, t};
Clear[a, b, c, d];
With[{a = a[t, y], b = b[t, y], c = c[t, y], d = d[t, y]},
eq = {I*a + D[b, y] + I*c == 0,
D[a, t] + I*f[y]*a + b*D[f[y], y] + I*d ==
1/m*(D[a, {y, 2}] - 2*a),
D[b, t] + I*f[y]*b + D[d, y] == 1/m*(D[b, {y, 2}] - 2*b),
D[c, t] + I*f[y]*c + I*d == 1/m*(D[c, {y, 2}] - 2*c)};
ic = {a == a0[y], b == b0[y], c == c0[y], d == d0[y]} /. t -> 0;
bc = {{a == 0, b == 0, c == 0, d == 0} /.
y -> -1, {a == 0, b == 0, c == 0, d == 0} /. y -> 1};];
ptoafunc = pdetoae[{a, b, c, d}[t, y], grid /@ {t, y}, difforder];
del = #[[2 ;; -2]] &;
ae = Map[del, Most /@ ptoafunc[eq], {2}];
aeic = Map[del, ptoafunc[ic]];
aebc = ptoafunc@bc;
var = Outer[#[#2, #3] &, {a, b, c, d}, grid[t], grid[y]];
{barray, marray} =
CoefficientArrays[Flatten[{ae, aeic, aebc}], Flatten[var]];
sollst = LinearSolve[marray, -barray];
solmatlst = ArrayReshape[sollst, var // Dimensions];
solfunclst = ListInterpolation[#, grid /@ {t, y}] & /@ solmatlst;
Update:
pdetoae can solve this system, however, the solution obtained is totally noisy, like the following figure. Also I found if using ic = {a == a0[y], b == b0[y], c == c0[y], d == d0[y]} /. t -> Tend, that is, using Tend as the initial time, the solution seems reasonable.
I am quite sure that the problem stems from some mistakes in the removal of boundary difference equations for implementation of BCs/ICs instead of using higher difference order. According to this answer, I understood that: the b.c.s of the problem are Dirichlet type so there's no need to use pdetoae on them, as did in the code.
Some specific questions:
In
ae = Map[del, Most /@ ptoafunc[eq], {2}], why we should applydelto each element on the 2nd level in listfptoafunc[feq]with its last element removed (byMost[list])?How to store the solution, which can then be used to solve for another equation? For example,
NDSolve[{D[u[t, y], t] - 1/m*D[u[t, y], {y, 2}] == -2*Real[b[t, y]*Conjugate[D[a[t, y], y]] + I*Conjugate[c[t, y]]*a[t, y]], u[0, y] == 0, u[t, -1] == 0, u[t, 1] == 0}, u, {t, 0, Tend}, {y, -1, 1}]
in which the solution a[t,y], b[t, y] and c[t, y] should be used to compute the RHS of the diffusion equation.
Could anybody help me? Thank you for any suggestion.
pdetoaeis merely auxiliary function for the implementation of FDM, you must first understand the basic of FDM before usingpdetoae. 2. "I didn't want to use Dumpsave since it normally is used to save the full info. " What do you mean by "full info"? $\endgroup$Most /@ ptoafunc[eq]although I knew the basic FDM. 2. By ''full info.'' I meant all data of anInterpolationFunction, including mesh, all the derivative values, etc. Here, only the time evolution of the original functions and the first space-derivatives are needed. If my understanding is wrong, please kindly suggest. Thanks in advance. $\endgroup$pdetoaefirst. 2. Your guess forInterpolatingFunctionis wrong. Just checkInterpolation[{1, 3, 6, 4, 5}] // InputForm. Also, check this post. $\endgroup$LinearSolve. " No, even if it's an explicit scheme, one can still write down every equation and solve the whole system all in once withLinearSolve. Yes, the answer for your question is replaceMostwithRest, now tell me why. $\endgroup$InterpolatingFunctionisn't as expensive as you've imagined. Store it withDumpSaveis the proper way to go. $\endgroup$