Change of variables in 2-nd order PDE

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?

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

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_] :=
NDSolveFiniteDifferenceDerivative[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 = DeveloperToPackedArray[L1, Real]; L2 =
DeveloperToPackedArray[L2, Real]; L3 =
DeveloperToPackedArray[L3, Real]; m1 =
DeveloperToPackedArray[m1, Real]; m2 =
DeveloperToPackedArray[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}]


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_] :=
NDSolveFiniteDifferenceDerivative[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 = DeveloperToPackedArray[L1, Real]; L2 =
DeveloperToPackedArray[L2, Real]; L3 =
DeveloperToPackedArray[L3, Real]; m1 =
DeveloperToPackedArray[m1, Real]; m2 =
DeveloperToPackedArray[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}]