NonStandard finite difference for BMBB equation

Question

I am trying to implement the following scheme mentioned in the paper

"NUMERICAL SOLUTIONS OF BENJAMIN-BONA-MAHONY-BURGERS EQUATION VIA NONSTANDARD FINITE DIFFERENCE SCHEME " What is mentioned in the paper for the BMMB Equation : http://math-frac.org/Journals/EJMAA/Vol6(2)_July_2018/Vol6(2)_Papers/21_EJMAA_Vol6(2)_July_2018_pp_237-245.pdf

Using the following scheme

I have tried this code in Mathematica

How to discretize a nonlinear PDE fast? nonlinear-pde-fast/28011 Code Edited following the post

NN = 8 ;
M = 8 ;
a = -10 ; 
b = 10 ;
h = (b - a)/NN;
T = 0.5 ;
k = T/M ; 
\[Phi][x_] = (E^(Sqrt[2] h/3) - 2 + E^(-Sqrt[2] h/3))/(2/9);
\[Psi][y_] = Sinh[y];
(*Defining the Grid points*)
Table[Subscript[x,i] = -10 + i h , {i , 0, M}];
Table[Subscript[t,j] = 0 + j k , {j , 0 , NN } ] ;
(*Defining the Initial Conditions*)
For[i = 0, i <= M, i++ ,Subscript[w , i , 0 ] = E^(-Subscript[x, i]*Subscript[x, i])];
(*Defining the Boundary Conditions*)
For [j = 0 , j <= NN , j++, Subscript[w, 0, j]  = 0];
For[j = 0 , j <= NN , j++ ,Subscript[w, M, j]  = 0];
(Defining the nonlinear equations due to discritization)
For[i = 1 , i <= NN, 
  i++ , {For [j = 1, j <= M - 1 , j++, 
    f[i, j] = 
     Subscript[w, i + 1, 
        j + 1](1/(2[Psi][k]([Phi][h])^2) + 1/(
         2 ([Phi][h])^2) - 
         1/(4[Phi][h])(1 + (
            Subscript[w, i, j + 1] + Subscript[w, i, j])/2))
      - ((Subscript[w, i, j + 1] - Subscript[w, i, j - 1])/(
        2* [Psi][k]) - (-2 Subscript[w, i, j + 1] + 
         Subscript[w, i - 1, j + 1] - Subscript[w, i + 1, j - 1] + 
         2Subscript[w, i, j - 1] - Subscript[w, i - 1, j - 1])/(
        2* [Psi][k]([Phi][h])^2)
        - (-2Subscript[w, i, j + 1] + Subscript[w, i - 1, j + 1] + 
         Subscript[w, i + 1, j - 1] - 2Subscript[w, i, j - 1] + 
         Subscript[w, i - 1, j - 1])/(2([Phi][h])^2
        )
        + (1 + (Subscript[w, i, j + 1] + Subscript[w, i, j])/
           2)((-Subscript[w, i - 1, j + 1] + Subscript[w, i + 1, j] -
            Subscript[w, i - 1, j])/(4([Phi][h]))))]}];
Sys = Flatten[Table[f[i, j], {i, M - 1}, {j, NN }]]//FullSimplify;
Vec = Flatten[Table[Subscript[w, i, j], {i, M - 1}, {j, NN}]];
Sol = FindRoot[Sys, {#, 1} & /@ Vec]

I get an error with FindRoot
How can i fix my code to show the result For N=M=8 ?

0.I am trying to implement the method mentioned in the paper .Discretization of time and space using central finite difference .1. Yes .I get an Error with FindRoot .How can i fix the error ? I think If i use N =8 , M =8 . I get 8 points out of the domain So i need to replace them with 0 to fix the problem . How can i replace all 8 values in the Sys ? — Mahmoud Hassan, Aug 13 '21 at 09:31
It's the paper that's unclear at this point. As you've noticed, with the scheme in the paper we still miss NN-1 equation, because the scheme uses non-standard central difference formula in $t$ direction, too. In traditional finite difference method, the missing equation can be supplied by one-sided difference formula, perhaps the author of the paper has done something similar, but this doesn't seem to be explained in the paper. — xzczd, Aug 13 '21 at 12:01

Alex Trounev · Accepted Answer · 2021-08-18T18:22:21.153

As it mentioned by xzczd we need to add some initial condition to this model. In a case of Example 2 we can use for example $u_t=0$ at $t=0$, then we have

u0[x_] := Exp[-x^2]; u1[x_] := 0;ue[x_,t_]:=0;
a = -10;
b = 10; NN = 200;
M = 50; alpha = 1; beta = 1; c = 1/10; v = 1 + c; kap = Sqrt[c/(4 v)];
h = (b - a)/NN; xcol = Table[a + h i, {i, 0, NN}];
T = 1/2; dt = T/M; tcol = Table[dt j, {j, 0, M}];
\[Phi][x_] := (E^(Sqrt[2] x/3) - 2 + E^(-Sqrt[2] x/3))/(2/9);
\[Psi][y_] := Sinh[y];
eqn = Table[(u[n, m + 1] - u[n, m - 1])/(2 \[Psi][dt]) - (u[n + 1, m + 1] - 
2 u[n, m + 1] + u[n - 1, m + 1] - u[n + 1, m - 1] + 
2 u[n, m - 1] - u[n - 1, m - 1])/(2 \[Psi][dt] \[Phi][h]) - 

alpha (u[n + 1, m + 1] - 2 u[n, m + 1] + u[n - 1, m + 1] + 
     u[n + 1, m - 1] - 2 u[n, m - 1] + 
     u[n - 1, m - 1])/(2 [Phi][h]) + (beta + (u[n, m + 1] + u[n, m])/
     2) (u[n + 1, m + 1] - u[n - 1, m + 1] + u[n + 1, m] - 
     u[n - 1, m])/(4 [Phi][h]^.5), {n, NN - 1}, {m, M - 1}];
    sys = Join[Flatten[Table[eqn[[i, j]] == 0, {i, NN - 1}, {j, M - 1}]], 
       Flatten[
        Table[{u[0, j] == ue[a, tcol[[j]]], 
          u[NN, j] == ue[b, tcol[[j]]]}, {j, 1, M}]], 
       Table[u[i, 0] == u0[xcol[[i + 1]]], {i, 0, NN}], 
       Table[u[i, 1] == u[i, 0] + dt u1[xcol[[i + 1]]], {i, 1, NN - 1}]];
    guess = Flatten[Table[{u[i, j], 1/10}, {i, 0, NN}, {j, 0, M}], 1];
sol = FindRoot[sys, guess];

Visualization

lst = Table[{xcol[[i + 1]], tcol[[j + 1]], u[i, j] /. sol}, {i, 0, 
    NN}, {j, 0, M}];
ListPlot3D[Flatten[lst, 1], PlotRange -> All, Mesh -> None, 
 ColorFunction -> "Rainbow", PlotTheme -> "Marketing"]

I have tested this numerical algorithm in a case of Example 1 and 3 as well. It is not so perfect as it mentioned in the paper. Also there are several typos in the equation (5) and (6) in the paper. I try to fix all of them. Code for Example 1:

u0[x_] := 3 c Sech[kap x]^2; ue[x_, t_] := 3 c Sech[kap (x - v t)]^2; 
u1[x_] := 3 c kap v Sech[kap x] Tanh[kap x];
a = -40;
b = 60; NN = 200;
M = 20; alpha = 0; beta = 1; c = 1/10; v = 1 + c; kap = Sqrt[c/(4 v)];
h = (b - a)/NN; xcol = Table[a + h i, {i, 0, NN}];
T = 1; dt = T/M; tcol = Table[dt j, {j, 0, M}];
\[Phi][x_] := (E^(Sqrt[2] x/3) - 2 + E^(-Sqrt[2] x/3))/(2/9);
\[Psi][y_] := Sinh[y];
eqn = Table[(u[n, m + 1] - 
       u[n, m - 1])/(2 [Psi][dt]) - (u[n + 1, m + 1] - 
       2 u[n, m + 1] + u[n - 1, m + 1] - u[n + 1, m - 1] + 
       2 u[n, m - 1] - u[n - 1, m - 1])/(2 [Psi][dt] [Phi][h]) - 
    alpha (u[n + 1, m + 1] - 2 u[n, m + 1] + u[n - 1, m + 1] + 
        u[n + 1, m - 1] - 2 u[n, m - 1] + 
        u[n - 1, m - 1])/(2 [Phi][
         h]) + (beta + (u[n, m + 1] + u[n, m])/2) (u[n + 1, m + 1] - 
        u[n - 1, m + 1] + u[n + 1, m] - 
        u[n - 1, m])/(4 [Phi][h]^.5), {n, NN - 1}, {m, M - 1}];
sys = Join[Flatten[Table[eqn[[i, j]] == 0, {i, NN - 1}, {j, M - 1}]], 
   Flatten[
    Table[{u[0, j] == ue[a, tcol[[j]]], 
      u[NN, j] == ue[b, tcol[[j]]]}, {j, 1, M}]], 
   Table[u[i, 0] == u0[xcol[[i + 1]]], {i, 0, NN}], 
   Table[u[i, 1] == u[i, 0] + dt u1[xcol[[i + 1]]], {i, 1, NN - 1}]];
guess = Flatten[Table[{u[i, j], 1/10}, {i, 0, NN}, {j, 0, M}], 1];
sol = FindRoot[sys, guess];
lst = Table[{xcol[[i + 1]], tcol[[j + 1]], u[i, j] /. sol}, {i, 0, 
    NN}, {j, 0, M}];

Visualization of error

error = Table[{xcol[[i + 1]], tcol[[j + 1]], 
    Abs[ue[xcol[[i + 1]], tcol[[j + 1]]] - u[i, j] /. sol]}, {i, 0, 
    NN}, {j, 0, M}];
ListPointPlot3D[error, PlotTheme -> "Marketing", PlotRange -> All]

We can compare this result with NDSolve using code @xzczd

With[{u = u[t, x]}, 
 eq = D[u, t] - D[u, x, x, t] - 
    alpha D[u, x, x] + (beta + u) D[u, x] == 0; 
 ic = u == u0[x] /. t -> 0; 
 bc = u == u0[x] /. {{x -> a}, {x -> b}}]; tst2 = 
 Interpolation[Flatten[lst, 1]]; nsol = 
 NDSolveValue[{eq, ic, bc}, u, {t, 0, 1}, {x, a, b}];
With[{t = 1}, 
 Plot[{tst2[x, t], nsol[t, x]}, {x, a, b}, PlotRange -> All]]
With[{t = 1}, 
 Plot[{Abs[tst2[x, t] - ue[x, t]], Abs[nsol[t, x] - ue[x, t]]}, {x, a,
    b}, PlotRange -> All]]

The result is inconsistent with that of NDSolve: With[{u = u[t, x]}, eq = D[u, t] - D[u, x, x, t] - alpha D[u, x, x] + (beta + u) D[u, x] == 0; ic = u == Exp[-x^2] /. t -> 0; bc = u == 0 /. {{x -> -10}, {x -> 10}}]; tst2 = Interpolation[Flatten[lst, 1]]; nsol = NDSolveValue[{eq, ic, bc}, u, {t, 0, 0.5}, {x, -10, 10}]; With[{t = 0.5}, Plot[{tst2[x, t], nsol[t, x]}, {x, -10, 10}, PlotRange -> All]] . Not sure what's wrong… — xzczd, Aug 14 '21 at 08:17
@MahmoudHassan See update to my answer with code for Example 1. — Alex Trounev, Aug 18 '21 at 16:41
@xzczd There are several typos in the paper in equation (6), I have fixed some of them and now the results are consistent. — Alex Trounev, Aug 18 '21 at 18:17
@AlexTrounev I didn't notice . What are the typos in the paper ? — Mahmoud Hassan, Aug 21 '21 at 07:56
Why ? u1[x_] := 3 c kap v Sech[kap x] Tanh[kap x]; Instead of u1[x_] := 6 c kap v Sech[kap x]^2 Tanh[kap x]; — Mahmoud Hassan, Aug 21 '21 at 07:56
And Why did you decreased the exponent for [Phi][h]^2 you wrote [Phi][h] And instead of [Phi][h] you wrote [Phi][h]^0.5 ? — Mahmoud Hassan, Aug 21 '21 at 08:04
@MahmoudHassan They suppose first in the paper that $\phi (h)=h+O(h^2)$ while then they use $\phi (h)=h^2+O(h^4)$. We can change code for u[n, m] as I did or change $\phi (h)$. Also I have used 3 instead of 6 in u1[x] to minimize error. — Alex Trounev, Aug 21 '21 at 09:54

NonStandard finite difference for BMBB equation

1 Answers1