3
$\begingroup$

I need to solve the following equation:

$-(1+x) \partial_t ^2 \psi + (1-2x^2)\partial_x\partial_t \psi +(1-x)x^2\partial_x^2 \psi -2x \partial_t\psi+x(2-3x)\partial_x\psi=0$

for the function $\psi(t,x)$, for $x \in [0,1]$, with initial conditions $\psi(0,x)=f(x), \dot{\psi}(0,x)=(1-x)f'(x)$, where $f(x)=e^{-218(\frac{-2}{3}-ln(1-x))^2}$. The idea is to use the change of variables $\pi = \partial_t \psi$ to reduce it to a system that is 1-st order in time, $\partial_t u=\hat{L}u$, where $u=\begin{bmatrix} \psi \\ \pi \end{bmatrix}, \hat{L}= \begin{bmatrix} 0 && 1 \\ \hat{L_1} && \hat{L_2} \end{bmatrix}$. The system will then be solved with a pseudo-spectral method for spatial discretization and a Runge-Kutta method for time integration, in the spirit of this question, but I'm having trouble with defining the new equation and the matrix $\hat{L}$ under the change of variables, in order to use it in the latter methods, can someone offer a proposal?

Improve this question
$\endgroup$
10
  • 2
    $\begingroup$ Have you seen DChange[]? It might work on this. $\endgroup$ Commented Jan 19, 2022 at 17:31
  • $\begingroup$ What are interval of integration on x and bc? $\endgroup$ Commented Jan 19, 2022 at 17:47
  • $\begingroup$ @AlexTrounev Sorry, I edited the question to include them $\endgroup$ Commented Jan 19, 2022 at 18:00
  • $\begingroup$ @JBuck Boundary conditions are not defined. Should we integrate without bc with ic only? $\endgroup$ Commented Jan 19, 2022 at 23:44
  • $\begingroup$ Why not directly use NDSolve? Is it for didactic purposes? $\endgroup$ Commented Jan 20, 2022 at 0:42

1 Answer 1

Reset to default
4
$\begingroup$

To solve this problem we can use method of line with handmade discreditation on x as follows:

Clear["Global`*"] Needs["Developer`"] b = -2/3; a = 218; x0 = 0; x1 = 1; u = Exp[-a*(b + t - Log[1 - x])^2]; f[y_] := u /. {t -> 0, x -> y}; f1[y_] := D[u, x] /. {t -> 0, x -> y}; eq = -(1 + x) D[u, t, t] + (1 - 2 x^2)* D[u, x, t] + (1 - x) x^2 D[u, x, x] - 2 x D[u, t] + x (2 - 3 x) D[u, x] // Simplify der[x0_, k_] := NDSolve`FiniteDifferenceDerivative[Derivative[k], x0, "DifferenceOrder" -> "Pseudospectral"]@"DifferentiationMatrix"; Nx = 2^8; II = IdentityMatrix[Nx]; n = Nx - 1; theta = N[Table[i*Pi/n, {i, 0, n}]]; X = N[Chop[((x0 + x1)/2) + ((x0 - x1)/2)*Cos[theta]]]; D1 = der[X, 1]; D2 = der[X, 2]; L1 = (II - 2 DiagonalMatrix[X]^2) . D1; L2 = ((II - DiagonalMatrix[X]) . DiagonalMatrix[X]^2) . D2; L3 = (DiagonalMatrix[X] . (2 II - 3 DiagonalMatrix[X])) . D1; m1 = II + DiagonalMatrix[X]; m2 = 2 DiagonalMatrix[X]; L1 = Developer`ToPackedArray[L1, Real]; L2 = Developer`ToPackedArray[L2, Real]; L3 = Developer`ToPackedArray[L3, Real]; m1 = Developer`ToPackedArray[m1, Real]; m2 = Developer`ToPackedArray[m2, Real]; Psi = Table[Subscript[psi, i][t], {i, Nx}]; Psi1 = Table[Subscript[psi, i]'[t], {i, Nx}]; Pi0 = Table[Subscript[pi, i][t], {i, Nx}]; Pi1 = Table[Subscript[pi, i]'[t], {i, Nx}]; eq1 = -m1 . Pi1 + L1 . Pi0 + L2 . Psi - m2 . Pi0 + L3 . Psi; eq2 = -Psi1 + Pi0; F0 = Chop[f[X] /. Indeterminate -> 0] // Quiet; ic1 = Table[Subscript[psi, i][0] == F0[[i]], {i, Nx}]; F1 = (II - DiagonalMatrix[X]) . Chop[f1[X] /. Indeterminate -> 0] // Quiet; ic2 = Table[Subscript[pi, i][0] == F1[[i]], {i, Nx}]; var = Join[Table[Subscript[psi, i], {i, Nx}], Table[Subscript[pi, i], {i, Nx}]]; sol = NDSolve[Join[Table[eq1[[i]] == 0, {i, Nx}], Table[eq2[[i]] == 0, {i, Nx}], ic1, ic2], var, {t, 0, 1}];

Visualization of solution and absolute error

Table[{ListPlot[Transpose[{X, Psi /. sol[[1]] /. t -> tn}], PlotRange -> All, Joined -> True, PlotLabel -> tn], ListPlot[ Transpose[{X, Abs[(Psi /. sol[[1]] /. t -> tn) - (u /. {t -> tn, x -> X}) // Quiet]}], PlotRange -> All, Joined -> True, ScalingFunctions -> "Log"]}, {tn, 0.01, .91, .1}]

Figure 1

Update 1. We also can implement Runge-Kutta method from my answer here as follows

Clear["Global`*"] Needs["Developer`"] b = -2/3; a = 218; x0 = 0; x1 = 1; u = Exp[-a*(b + t - Log[1 - x])^2]; f[y_] := u /. {t -> 0, x -> y}; f1[y_] := D[u, x] /. {t -> 0, x -> y}; eq = -(1 + x) D[u, t, t] + (1 - 2 x^2)* D[u, x, t] + (1 - x) x^2 D[u, x, x] - 2 x D[u, t] + x (2 - 3 x) D[u, x] // Simplify der[x0_, k_] := NDSolve`FiniteDifferenceDerivative[Derivative[k], x0, "DifferenceOrder" -> "Pseudospectral"]@"DifferentiationMatrix"; Nx = 2^8; II = IdentityMatrix[Nx]; n = Nx - 1; theta = N[Table[i*Pi/n, {i, 0, n}]]; X = N[Chop[((x0 + x1)/2) + ((x0 - x1)/2)*Cos[theta]]]; D1 = der[X, 1]; D2 = der[X, 2]; L1 = (II - 2 DiagonalMatrix[X]^2) . D1; L2 = ((II - DiagonalMatrix[X]) . DiagonalMatrix[X]^2) . D2; L3 = (DiagonalMatrix[X] . (2 II - 3 DiagonalMatrix[X])) . D1; m1 = II + DiagonalMatrix[X]; m2 = 2 DiagonalMatrix[X]; L1 = Developer`ToPackedArray[L1, Real]; L2 = Developer`ToPackedArray[L2, Real]; L3 = Developer`ToPackedArray[L3, Real]; m1 = Developer`ToPackedArray[m1, Real]; m2 = Developer`ToPackedArray[m2, Real]; F0 = Chop[f[X] /. Indeterminate -> 0] // Quiet; F1 = (II - DiagonalMatrix[X]) . Chop[f1[X] /. Indeterminate -> 0] // Quiet; z = ConstantArray[0, 2 Nx]; ff[t_, z_] := Join[Drop[z, Nx], Inverse[m1] . (L1 . Drop[z, Nx] + L2 . Drop[z, -Nx] - m2 . Drop[z, Nx] + L3 . Drop[z, -Nx])]; rk2[ff_, h_][{t_, x_}] := Module[{k1, k2}, k1 = ff[t, x]; k2 = ff[t + h/2, x + h k1/2]; {t + h, x + h k2}] tf = 1; dt = 1/1000; sol = NestList[rk2[ff, dt], {0, Join[F0, F1]}, Round[tf/dt]];

Visualization of numerical solution (red points) and exact solution (blue lines) at different times (shown above)

Table[Show[ ListLinePlot[Transpose[{X, (u /. {t -> i dt, x -> X}) // Quiet}], PlotStyle -> Blue, PlotRange -> All, PlotLabel -> N[i dt]], ListPlot[Transpose[{X, Drop[sol[[i + 1, 2]], -Nx]}], PlotStyle -> Red, PlotMarkers -> {Automatic, Smaller}]], {i, 1, Round[tf/dt], 100}]

Figure 2

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.