I want to solve two coupled PDEs with pdetoode function:
feq = {D[ϕ, t] + U*D[ϕ, y] - D[ϕ, {y, 2}] - (s0 - 1/20*y^2)*ϕ == 0}
and
beq = {D[a, t] + U*D[a, y] + D[a, {y, 2}] + (s0 - 1/20*y^2)*a == 0}
To this end, I encapsulated them in the two solvers, which are very similar
(function pdetoode not included here, please find it from the above link):
feqSolver[s_, T_, g_] := Module[{s0 = s, Tend = T, feqIC = g},
With[{ϕ = ϕ[t, y], U = y},
feq = {D[ϕ, t] + U*D[ϕ, y] -
D[ϕ, {y, 2}] - (s0 - 1/20*y^2)*ϕ == 0};
fic = {ϕ == feqIC[y]} /. t -> 0;
fbc = {{ϕ == 0} /. y -> Ly, {ϕ == 0} /. y -> Uy};];
fptoofunc = pdetoode[ϕ[t, y], t, grid, difforder];
fdel = #[[2 ;; -2]] &;
fode = fdel@fptoofunc@feq[[1]];
fodeic = fptoofunc@fic;
fodebc = fptoofunc@fbc;
fsollst =
NDSolveValue[{fodebc, fodeic, fode},
Map[ϕ, grid], {t, 0, Tend}];
fsol = ListInterpolation[
Developer`ToPackedArray@#["ValuesOnGrid"] & /@ fsollst //
Transpose, {Flatten@fsollst[[1]]["Grid"], grid}];
phisol[t_, y_] = ϕ[t, y] /. ϕ -> fsol[[1]];]
beqSolver[s_, T_, aT_] :=
Module[{s0 = s, Tend = T, beqIC = aT},
With[{a = a[t, y], U = y},
beq = {D[a, t] + U*D[a, y] + D[a, {y, 2}] + (s0 - 1/20*y^2)*a == 0};
bic = {a == beqIC[y]} /. t -> Tend;
bbc = {{a == 0} /. y -> Ly, {a == 0} /. y -> Uy};];
bptoofunc = pdetoode[a[t, y], t, grid, difforder];
bdel = #[[2 ;; -2]] &;
bode = bdel@bptoofunc@beq[[1]];
bodeic = bptoofunc@bic;
bodebc = bptoofunc@bbc;
bsollst =
NDSolveValue[{bodebc, bodeic, bode}, Map[a, grid], {t, Tend, 0}];
bsol = ListInterpolation[
Developer`ToPackedArray@#["ValuesOnGrid"] & /@ bsollst //
Transpose, {Flatten@bsollst[[1]]["Grid"], grid}];
asol[t_, y_] = a[t, y] /. a -> bsol[[1]];]
Spatial domain and initial condition (IC):
Clear[g0, Ly, Uy]
Ly = -20; Uy = 20;
(*Giving a random IC for feq*)
SetAttributes[g0, Listable];
SeedRandom[1];
g0[y_?NumericQ] =
BSplineFunction[Join[{0.}, RandomReal[{-1, 1}, 39], {0.}],
SplineClosed -> False][(y + 20)/40];
Define two functions: J is a function of interest, aT is a function used for constructing the IC for beq
Clear[J, aT];
J[g_, phiT_, y_] :=
NIntegrate[Dot[g[y], g[y]], {y, Ly, Uy}]/
NIntegrate[Dot[phiT[y], phiT[y]], {y, Ly, Uy}];
(*define how to generate IC for `beq` from the soln of `feq`*)
aT[g_, phiT_,
y_] := (-2*phiT[y]*
NIntegrate[Dot[g[y], g[y]], {y, Ly, Uy}])/(NIntegrate[
Dot[phiT[y], phiT[y]], {y, Ly, Uy}])^2;
Solve the two equation:
s0 = 0.4; Tend = 10;
g[y_?NumericQ] = g0[y];
phisol[t, y] = feqSolver[s0, Tend, g];
phiT[y] = phisol[Tend, y];
Jtest[[1]] = J[g, phiT, y]; (*evaluate a function of interest and save it as the first element of a list Jtest*)
aT[y] = aT[g, phiT, y];
asol[t, y] = beqSolver[s0, Tend, aT];
a0[y] = asol[0, y];
When run this code, MMA mainly gives this kind of error:
Part::take: "Cannot take positions 2 through -2 in pdetoode[ϕ[t,y],t,grid,difforder]" $RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>
By running the code line by line, I found that the first error results from phisol[t, y] = feqSolver[s0, Tend, g];
In the above code, I have used Jtest[[1]] since in my real problem I want to evaluate and save Jtest[[iter]] continuously (in a list) in a loop with an iterator iter. Then Jtest[[iter]] will be used for stopping the loop as follows:
Jtest[[iter]] = J[g, phiT, y]; AppendTo[Jtest, Jtest[[iter]]];
If[Abs[Jtest[[iter]] - Jtest[[iter - 1]]] < acc, Break[]];
**My problems are (1) should I define an empty list in the beginning, i.e. something like Jtest={}?; and (2) how do I fix these problems?
I have searched on the SE site for a similar problem but failed to solve it so far. And I have tried my best to learn the fundamental syntax of Mathematica for several months but I'm still in trouble. Please help and any suggestions and comments are very welcome. Thank you in advance.
J[[iter]] = J[g, phiT, y];doesn't make sense at all. Please consider asking a question about this, rather than immediately trying to deal with such a big problem. Also, recently I've written a relevant Chinese tutorial, please consider reading it: https://note.youdao.com/share/?id=f2cdd9e72b66c73f372296800de7f7f5&type=note#/ – xzczd Dec 20 '18 at 10:21beq = {D[a, t] + U*D[a, y] + D[a, {y, 2}] + (s0 - 1/20*y^2)*a == 0}? I ask about the sign+ D[a, {y, 2}]. – Alex Trounev Dec 20 '18 at 16:20pdetoode/pdetoaeis just auxiliary tool for implementation of FDM, so it can of course handle PDE system, see e.g. https://mathematica.stackexchange.com/a/184285/1871 https://mathematica.stackexchange.com/q/278430/1871 https://mathematica.stackexchange.com/q/183745/1871 And you'll find more by searchingpdetoode/pdetoaein this site. – xzczd Oct 31 '23 at 01:41