1

I have been trying to solve a partial differential equation known as KdV (Korteweg-de Vries) equation using NDSolve.
KdV equation is written as:

$\frac{\partial u(x,t)}{\partial t}+\frac{\partial^3 u(x,t)}{\partial x^3}=-6u(x,t)\frac{\partial u(x,t)}{\partial x}$

One of the analytical solution of this equation is a Jacobi elliptical function
$u_{A} = \text{A JacobiCN}[B*(x- c t),m]^2$
where A, B, c and m are known parameters.

So the initial condition is
$u_0 = \text{A JacobiCN}[B*x,m]^2 $.

The boundary condition is: $u(t,-a) = u(t,a)$.

Code for analytical solution is

a = 30; m = 0.9; A = 0.9;
B = Sqrt[1.5*A] / Sqrt[2*m + m^2]
c = 6*(-A + A*m + A*m^2)/(m*(2 + m))
uA[x_,t_] := A*JacobiCN[B*(x - c*t), m]^2;
Plot[uA[x,0.5], {x, -a, a}, LabelStyle ->  Directive[14, Bold,   Black], AxesLabel -> {Style["x", 14, Bold, Black], Style["uA", 14, Bold, Black]}]

Plot of analytical solution for t = 0.5

enter image description here

When I solve the partial differential equation using NDSolve, the numerical instabilities arise and the numerical solution doesn't match the analytical solution.

eq = {D[u[x, t], {t, 1}] + D[u[x, t], {x, 3}] == -6 u[x, t] D[u[x, t], {x, 1}], u[x, 0] == A*JacobiCN[B*x, m]^2, u[-a, t] == u[a, t]};
sol = Flatten@NDSolve[eq, u, {t, 0, 1}, {x, -a, a}]
Plot[Evaluate[u[x, 0.5] /. sol], {x, -a, a}, LabelStyle ->  Directive[14, Bold, Black], AxesLabel -> {Style["x", 14, Bold, Black], Style["u", 14, Bold, Black]}]

enter image description here

I am also getting the warning message NDSolve::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x. I have tried different options and methods given in NDSolve documentation and also tried solutions given in https://mathematica.stackexchange.com/questions/9277/i-failed-to-solve-a-set-of-one-dimension-fluid-mechanics-pdes-with-ndsolve/9291#9291 but nothing seems to be working. Following method removed the warning message but there was no improvement in the numerical solution.

nxy = 500;
sol = Flatten@NDSolve[eq, u, {t, 0, 1}, {x, -a, a}, Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> nxy, "MinPoints" -> nxy, 
       "DifferenceOrder" -> "Pseudospectral"}, Method -> "Adams"}, MaxSteps -> \[Infinity]]
bt89
  • 11
  • 1
  • 2
    The uA isn't an analytic solution of eq, just try Subtract @@@ eq /. u -> uA /. x -> a/2 /. t -> 5. Please double check your equation system. – xzczd Oct 20 '20 at 06:29
  • It is an exact solution. If I substitute the analytical solution into the partial differential equation and simplify it, I will get a common factor. When I write JacobiDN and JacobiSN in terms of JacobiCN, this common factor goes to zero. – bt89 Oct 21 '20 at 06:35
  • As shown in my last comment, numeric calculation indicates the uA doesn't satisfy the system. What do you mean by "If I substitute the analytical solution into the partial differential equation and simplify it, I will get a common factor. When I write JacobiDN and JacobiSN in terms of JacobiCN, this common factor goes to zero"? Can you show the code for this verification? – xzczd Oct 21 '20 at 06:49
  • u[x_, t_] := A*JacobiCN[B*(x - C1*t), m]^2; eq = Simplify[ D[u[x, t], {t, 1}] + D[u[x, t], {x, 3}] + 6 u[x, t]*D[u[x, t], {x, 1}]];f1 = (1.9016506529947896 - 4.5800318543959JacobiCN[1.0564725849971053 t - 0.7191949522280763x, 0.9]^2 + 2.678381201401112JacobiDN[ 1.0564725849971053 t - 0.7191949522280763x, 0.9]^2 - 2.4105430812610016JacobiSN[1.0564725849971053 t - 0.7191949522280763x, 0.9]^2);` – bt89 Oct 21 '20 at 07:23
  • Simplify[f1 /. {JacobiSN[ 1.0564725849971053 t - 0.7191949522280763x, 0.9]^2 -> 1 - JacobiCN[1.0564725849971053t - 0.7191949522280763 x, 0.9`]^2, JacobiDN[1.0564725849971053t - 0.7191949522280763 x, 0.9`]^2 -> ((1 - m^2) + m^2 JacobiCN[1.0564725849971053t - 0.7191949522280763 x, 0.9`]^2)}]` Sorry I was having difficulty in writing the code in a single comment so I had to break it. f1 is the common factor you get after simplifying eq. – bt89 Oct 21 '20 at 07:27
  • Please add the code to your question by clicking the edit button at the bottom left corner of the question. Then, 2nd replacing rule is wrong, try: JacobiDN[1.0564725849971053 t - 0.7191949522280763` x, 0.9]^2 -> ((1 - m^2) + m^2 JacobiCN[1.0564725849971053 t - 0.7191949522280763` x, 0.9]^2) /. t -> 1/2 /. x -> 2 – xzczd Oct 21 '20 at 07:42

0 Answers0