I am using differential quadrature method (DQM) to solve the following initial-boundary value problem:
$$\frac{\partial^2 w}{\partial t^2} + \mu \frac{\partial w}{\partial t} + C_0 \frac{\partial^4 w}{\partial x^4}=0$$ $$w(0,t)=X_b(t)$$ $$\frac{\partial w}{\partial x}\bigg|_{x=0}=\frac{\partial^2 w}{\partial x^2}\bigg|_{x=1}=\frac{\partial^3 w}{\partial x^3}\bigg|_{x=1}=0$$
Here $X_b(t)$ (Xb[t]
in the code below) is the input in the system. For harmonic function
$$X_b(t)=2G_1 \cos(2\Omega t)$$
as input, NDSolve
works perfectly. For other inputs also, simulation goes on properly. But for the input
$$X_b(t)=2G \cos(2 \Omega t) (w(1,t)+\alpha \; w(1,t)^3)$$
High frequency oscillations grow and the simulation results become more and more inaccurate, as number of grid points (Np
in the code below) increases. If tmax
is large or Np > 10
, then the simulation fails with warning
NDSolve::ndsz: singularity or stiff system suspected.
I have solved this problem with another method, there is no high frequency oscillations present.
The following is my trial. The PDE has been discretized to Np - 1
ODEs with DQM.
Np = 10; G1 = 0.05; Ω = 1; μ = 0.05; tmax = 10; a = 30;
ii = Range[1, Np]; X = 0.5 (1 - Cos[(ii - 1)/(Np - 1) π]);
Xi[xi_] := Cases[X, Except[xi]];
M[xi_] := Product[xi - Xi[xi][[l]], {l, 1, Np - 1}];
C1 = C3 = ConstantArray[0, {Np, Np, 4}];
Table[If[i != j, C1[[i, j, 1]] = M[X[[i]]]/((X[[i]] - X[[j]]) M[X[[j]]])];,
{i, 1, Np}, {j, 1, Np}];
Table[C1[[i, i, 1]] = -Total[Cases[C1[[i, All, 1]], Except[C1[[i, i, 1]]]]];,
{i, 1, Np}];
C3[[All, All, 1]] = C1[[All, All, 1]];
C3[[1, All, 1]] = 0;
C3[[All, All, 2]] = C1[[All, All, 1]].C3[[All, All, 1]];
C3[[Np, All, 2]] = 0;
C3[[All, All, 3]] = C1[[All, All, 1]].C3[[All, All, 2]];
C3[[Np, All, 3]] = 0;
C3[[All, All, 4]] = C1[[All, All, 1]].C3[[All, All, 3]];
C3r4 = N[C3[[All, All, 4]]];
C0 = 1.8751^-4;
K1M = C0 C3r4[[2 ;; Np, 1 ;; Np]];
Y1[t_] := Table[Subscript[x, i][t], {i, 2, Np}];
α = -0.001;
Input Xb[t]
is substituted in the system of equations through a column vector YV[t]
.
Xb[t] = 2 G1 Cos[2 Ω t] (Subscript[x, Np][t] + α Subscript[x, Np][t]^3);
YV[t] = Flatten[{Xb[t], Y1[t]}];
eqns = Thread[D[Y1[t], t, t] + μ D[Y1[t], t] + K1M.YV[t] == 0];
Lg = 1; bt = 1.8751/Lg; ξ = X[[2 ;; Np]];
y0 = -0.5 a (((Cos[bt*ξ] - Cosh[bt*ξ])-0.734*(Sin[bt*ξ] - Sinh[bt*ξ])));
X0 = -y0; X0d = 0 y0;
s = NDSolve[{eqns, Y1[0] == X0, Y1'[0] == X0d}, Y1[t], {t, 0, tmax},
Method -> {"StiffnessSwitching", Method->{"ExplicitRungeKutta", Automatic}},
MaxSteps -> Infinity,
AccuracyGoal -> 11, PrecisionGoal -> 11]; //AbsoluteTiming
plot1 = Plot[Evaluate[Subscript[x, Np][t] /. First@s], {t, 0, tmax},
PlotRange -> All]
Results obtained with version 11.1 for Np = 6
and 8
are given below. For Np = 8
, fluctuation in the output increases.
For Np = 12
, it gives warning
NDSolve::ndsz: At t == 7.129860412016928`, step size is effectively zero; singularity or stiff system suspected.
I have used different options of NDSolve
to deal with the stiff system, but still it is not working.
I thought that, if I use NDSolve
in smaller interval then it may be successful. So I made the code in which initial conditions (for ith iteration) are calculated based on the output from previous iteration (i - 1). But multiple NDSolve
simulation also didn't work.
I have tried different initial conditions also, but the warning remains. Can someone please help me to solve the issue? Thanks.
Plot[Evaluate[Subscript[x, Np][t] /. First@s], {t, 5.5, 6.5}, PlotRange -> {-70, 30}]
to see that a rapidly growing high frequency oscillation is superimposed on the slowly oscillating solution that you appear to desire. UsingWorkingPrecision -> 30
along withSetPrecision[{eqns, Y1[0] == X0, Y1'[0] == X0d}, 30]
does not help, which makes me think that the oscillation may be real. $\endgroup$WorkingPrecision -> 60
withRationalize[{eqns, Y1[0] == X0, Y1'[0] == X0d}, 0]
produces the same result! All this suggests that rapid, growing oscillations are real. $\endgroup$X[[0]]=0
, but have used all other points? Can you add explanation of the physical background of your code? It looks like nobody can understand what you try to solve. $\endgroup$