1
$\begingroup$

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 enter image description here

Using the following scheme enter image description here

enter image description here

enter image description here enter image description here 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] + 
         2*Subscript[w, i, j - 1] - Subscript[w, i - 1, j - 1])/(
        2* \[Psi][k]*(\[Phi][h])^2)
        - (-2*Subscript[w, i, j + 1] + Subscript[w, i - 1, j + 1] + 
         Subscript[w, i + 1, j - 1] - 2*Subscript[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 ?

$\endgroup$
3
  • 1
    $\begingroup$ 0. Why not NDSolve? 1. Have you read this post?: mathematica.stackexchange.com/q/10453/1871 2. m is undefined, please always pay attention to the color of the variable, the m is blue, which indicates it's "empty". 3. /FullSimplify is obviously wrong. 4. Think about what's wrong with the following: FindRoot[{x == 1, y == 2}, {{{x, 1}, {y, 2}}}] Check the document of Flatten and think about how to fix the sample with this function. 5. Check what's inside Guess[[1, 1]] and think about what's wrong. $\endgroup$
    – xzczd
    Aug 13, 2021 at 2:31
  • $\begingroup$ 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 ? $\endgroup$ Aug 13, 2021 at 9:31
  • 1
    $\begingroup$ 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. $\endgroup$
    – xzczd
    Aug 13, 2021 at 12:01

1 Answer 1

2
$\begingroup$

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"]

Figure 1

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]

Figure 2

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]]

Figure 3

$\endgroup$
11
  • $\begingroup$ 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… $\endgroup$
    – xzczd
    Aug 14, 2021 at 8:17
  • $\begingroup$ Thank you . For Example 1 What is the max error you got ? $\endgroup$ Aug 18, 2021 at 15:12
  • $\begingroup$ And What is the extra condition you added for example 1 ? $\endgroup$ Aug 18, 2021 at 15:49
  • $\begingroup$ @MahmoudHassan See update to my answer with code for Example 1. $\endgroup$ Aug 18, 2021 at 16:41
  • $\begingroup$ @xzczd There are several typos in the paper in equation (6), I have fixed some of them and now the results are consistent. $\endgroup$ Aug 18, 2021 at 18:17

Your Answer

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

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