Time & Space matlab discretization Finite Differences confusion

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);
end
end


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);
end


Function

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);

end
end
function du=d(u)

du = u(i+1)-u(i-1);
end
function du3 = d3(u)
du3 = u(i+2)-2*u(i+1)+2*u(i-1)-u(i-2);
end

u2 = u;
end


Image of what the plot should look like

• 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. Dec 6, 2021 at 11:29
• @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. Dec 6, 2021 at 15:58

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


or

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)

• 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 Dec 6, 2021 at 15:55
• I think that worked actually the only thing I changed was fu = 6*u.*du-du3 Dec 6, 2021 at 16:23
• I fixed the code thanks to you but I don't have enough reputation to upvote so I just wanted to say thanks. Dec 9, 2021 at 0:42