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?

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

1 Answer 1

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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