1
$\begingroup$

The problem I am trying to solve is the following :

$$(\frac{1}{c^2}\frac{\partial^2}{\partial t^2}-\frac{\partial^2}{\partial x^2}+1)\psi(x,t)=0$$

$$\psi (0,t)=\psi(L(t),t)=0$$

$$\psi(x,0)=\sqrt{\frac{1}{c\sqrt{1+\pi^2/4}}} sin(\pi \frac{x}{2})$$

$$\frac{d\psi (x,t)}{dt}\Bigr|_{t=0}=0$$

$$L(t)=2+sin(\omega t)$$

I have used the pdetoode method as was used here to solve this problem and here is my code (some codes might be redundant as I copied and modified the code from the link above) :

ydum = x/L[t];
eqn1 = 1/c^2 D[\[Psi][ydum, t], {t, 2}] - 
      D[\[Psi][ydum, t], {x, 2}] + (m^2 c^2)/h^2 \[Psi][ydum, t] == 
     0 /. \[Psi][ydum, t] -> \[Psi][y, t] /. x -> y L[t] // Expand

m = 1;
c = 1;
h = 1;
\[Omega] = 1;
L[t_] := 2 + Sin[\[Omega] t];
bc = {\[Psi][0, t] == 0, \[Psi][1, t] == 0};
ic = {\[Psi][y, 0] == 
    Sqrt[2 (m c)/(c Sqrt[m^2 c^2 + \[Pi]^2 h^2/L[0]^2] L[0])]
      Sin[ \[Pi] y], D[\[Psi][y, t], t] == 0 /. t -> 0};
points = 25;
xdifforder = 4;
{yL, yR} = domain = {0, 1};
grid = Array[# &, points, domain];
pto = pdetoode[\[Psi][y, t], t, grid, xdifforder];
remove = #[[2 ;; -2]] &;
{ode, odebc, odeic} = MapAt[remove, pto /@ {eqn1, bc, ic}, 1];
soll = NDSolveValue[{ode, odebc, odeic[[2]], odeic[[1]]}, \[Psi] /@ 
    grid, {t, 0, 50}];
sol1 = rebuild[#, grid, -1] &@soll
Manipulate[
 Plot[sol1[y, t] /. y -> y L[t], {y, 0, 1}, 
  PlotRange -> {-5, 10}], {t, 0, 50}]

enter image description here

The solution grows larger in time and this makes me anxious. Is this how it's supposed to be or is there something wrong with my code?

I have also solved this problem without the use of pdetoode (check out my previous question ) and I got different result : enter image description here

I greatly appreciate anyone who would point out my error in both methods.

$\endgroup$
3
  • 1
    $\begingroup$ The grid is too coarse, try e.g. points = 200; xdifforder = 2;. Also, it's better to observe a smaller time interval first e.g. {t, 0, 5}. $\endgroup$
    – xzczd
    Feb 19, 2022 at 13:36
  • $\begingroup$ Thanks! Now the solution doesn't get larger in time. However, when I ramp up the $\omega$ to 1000, NDSolveValue spits out an error : NDSolveValue::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. How do I fix this? (By the way, what does xdifforder mean in pdetoode?) $\endgroup$ Feb 20, 2022 at 5:18
  • 1
    $\begingroup$ 1. As to the icfail message, it's because you're currently discretize the system to a DAE system, to avoid it, you need to handle the b.c. in a more complicated way, see e.g. this post: mathematica.stackexchange.com/a/127411/1871 and this post mathematica.stackexchange.com/a/184285/1871 2. xdifforder sets the value of DifferenceOrder option of NDSolve`FiniteDifferenceDerivative, it's the order of difference formula for spatial discretization. $\endgroup$
    – xzczd
    Feb 20, 2022 at 14:31

0

Your Answer

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

Browse other questions tagged or ask your own question.