I have been trying to solve this equation and write the finite difference scheme in matlab for months, but I still am not successful.

Given the KdV Equation $$\tag{1}u_{t} -6uu_x+u_{xxx}=0$$ I have come up with a semi-discrete centered difference finite difference scheme which has been checked and confirmed by my professor to be.

$$\tag{2} \frac{\partial{u}}{\partial{t}}=6u(i)*\frac{u_{i+1}-u_{i-1}}{2\Delta x}-\frac{u_{i+2}-2u_{ui+1}+2u_{i-1}-u_{i-2}}{2\Delta x^3}$$

Given the solition equation $$\tag{3}u_1(x,t) = \frac{v}{2*cosh^2[1/2\sqrt{v}(x-vt-x0)]}$$

where $$v = 20 \ \ x0 = 0$$

I am also given $$\Delta x = 0.1 \ \ x \in [-10,10] \ \ \Delta t = 0.001 \ \ t_{span} = [0,2]$$ Periodic Boundary Conditions $$x(-10) = x(10)$$ Initial Conditions $$u(x,0) = u_1(x,0)$$ from eqn(3)

We are told to use a Runge_kutta scheme for the time integration $$\alpha_1 = \Delta t f(u^n)$$ $$\alpha_2 = \Delta t f(u^n+1/2\alpha_1)$$ $$\alpha_3 = \Delta t f(u^n+1/2\alpha_3)$$ $$\alpha_4 = \Delta t f(u^n+\alpha_4)$$ $$u^{n+1} = u^n + 1/6*(\alpha_1 + 2\alpha_2 + 2\alpha_3 + \alpha_4)$$ $$n = \text{time space} \ \ f = \text{spatial discretization function}$$

If somebody could please help me I don't know if it's my math that is wrong or my coding schemes. I've attached my code at the bottom if you want to take a look. I've only put one function that I've tried, but I have tried other things such as inverting the matrix before with no avail. After solving plot once at u(x,0) and then u(x,2)

``dx = 0.1; 
ga = (4*sqrt(6))/9; 
%dt = ga*dx^3; 
dt = 0.001; 
x1 = -10:dx:10; %x space 
t = 0:dt:2; %time 

v = 20; %solition thin g
x0 = 0; %x factor 

Nx = length(x1); 
Nt = length(t); 
x = linspace(-10,10,Nx+1); 
x(1) = []; 
%% Intiate Vectors
u = zeros(Nx,Nt); %u(x,t) space 
u1s = zeros(Nx,Nt); %solition solution 

%% Solition Solution 

for i = 1:Nx
    for j = 1:Nt
        u1s(i,j) = -v/(2*cosh(1/2*sqrt(v)*(x(i)-v*t(j)-x0))^2); 

Runge Kutta

u(:,1) = u0; 
ua = u; 
for j = 1:Nt
    a1 = dt*f5(u(:,j)); 
    a2 = dt*f5(u(:,j)+0.5*a1); 
    a3 = dt*f5(u(:,j)+1/2*a2); 
    a4 = dt*f5(u(:,j)+a3); 
%     u(1,j) = u(end,j);  %peridoic boundary 
    u(:,j+1) = u(:,j) + 1/6*(a1+2*a2+2*a3+a4); 


function u2 = f5(u)
%centered space-time discretization 

global dx dt Nx Nt

fu = 3*u.^2; 
c = dt/dx; 
h = dt/dx^3; 

for k = 0:dt:2
    for i = 3:Nx-2
        u = u + c/2*d(fu)-h/2*d3(fu); 
function du=d(u)
        du = u(i+1)-u(i-1); 
    function du3 = d3(u)
        du3 = u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2); 
u2 = u;

Image of what the plot should look like

enter image description here

  • $\begingroup$ KdV is only order 1 in time. There is another typo at the end of the discretization formula. Do you have any evidence that this numerical scheme with RK4 will work, and not blow up shortly after the start? The solution rapidly evolves ripples that get magnified to infinity. $\endgroup$ Dec 6, 2021 at 11:29
  • $\begingroup$ @LutzLehmann your right the discretization should be correct for $$u_t$$ I just edited the problem statement to be first order. As far as the Runge Kutta, that is what we are supposed to use in class. I was able to semi-solve it using a predictor corrector method but using RK as you said does create crazy ripples. $\endgroup$
    – bc_eng
    Dec 6, 2021 at 15:58

1 Answer 1


First you have to decide how to implement the periodic boundary conditions, as that influences the details of the state vector. For instance there is circshift that allows to implement the difference quotients as

um1 = circshift(u,1)
up1 = circshift(u,-1)

du = (up1-um1) ./ (2*dx)
d3u = (circshift(du,-1)+circshift(du,1)-2*du)/(dx*dx)

fu = 6*u.*du-d3u

Due to the periodicity, the value at x=10 duplicates the value at x=-10. Using this circular shift, the last value is provided by the first value where needed. Thus define

x1 = -10:dx:(10-dx); %x space 


x1 = -10:dx:10; %x space
x1 = x1(1:end-1); % last is duplicate of first 

The step-size has to satisfy $L\Delta t\le 0.5$ ($2.5$ is the stability boundary, but to get some accuracy usually this needs to be smaller), so with $L\simeq\frac1{\Delta x^3}$ this amounts to $Δt\le 0.5·Δx^3$, so from that side the chosen constants should work.

You should be able to use array operations for the reference solution

u1s = -v./(2*cosh((1/2*sqrt(v)).*(x'-v*t-x0)).^2); 

(assumption: x and t are row vectors, so that x' gets broadcast as identical columns and t as identical rows in x'-v*t)

  • $\begingroup$ Thank's for the tips I'll give it a try, but based on the stability conditions (I didn't post it as I thought it was unnecessary) we got $$\Delta T_{max} <= \frac{4\sqrt{6}}{9} \Delta x^3$$ which is about dt = 0.0011 $\endgroup$
    – bc_eng
    Dec 6, 2021 at 15:55
  • 1
    $\begingroup$ I think that worked actually the only thing I changed was fu = 6*u.*du-du3 $\endgroup$
    – bc_eng
    Dec 6, 2021 at 16:23
  • $\begingroup$ I fixed the code thanks to you but I don't have enough reputation to upvote so I just wanted to say thanks. $\endgroup$
    – bc_eng
    Dec 9, 2021 at 0:42

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.