Finite difference method for solving the Klein-Gordon equation in Matlab

Question

I am trying to numerically solve the Klein-Gordon equation that can be found here. To make sure I solved it correctly, I am comparing it with an analytical solution that can be found on the same link. I am using the finite difference method and Matlab. The initial spatial conditions are known, not the initial time conditions.

I start off by initializing the constants and the space-time coordinate system:

close all
clear 
clc
 
%% Constant parameters

A      = 2; 
B      = 3; 
lambda = 2; 
mu     = 3; 
a      = 4; 
b      = - (lambda^2 / a^2) + mu^2;  

%% Coordinate system

number_of_discrete_time_steps = 300; 
t = linspace(0, 2, number_of_discrete_time_steps); 
dt = t(2) - t(1); 

number_of_discrete_space_steps = 100; 
x = transpose( linspace(0, 1, number_of_discrete_space_steps) ); 
dx = x(2) - x(1); 

Next, I define and plot the analitical solution:

%% Analitical solution

Wa = cos(lambda * x) * ( A * cos(mu * t) + B * sin(mu * t) );

figure('Name', 'Analitical solution'); 
surface(t, x, Wa, 'edgecolor', 'none'); 
colormap(jet(256)); 
colorbar; 
xlabel('t'); 
ylabel('x'); 
title('Wa(x, t) - analitical solution');

The plot of the analytical solution is shown here.

In the end, I define the initial spatial conditions, execute the finite difference method algorithm and plot the solution:

%% Numerical solution

Wn = zeros(number_of_discrete_space_steps, number_of_discrete_time_steps);

Wn(1, :) = Wa(1, :); 
Wn(2, :) = Wa(2, :); 

for j = 2 : (number_of_discrete_time_steps - 1)  

    for i = 2 : (number_of_discrete_space_steps - 1) 

        Wn(i + 1, j) = dx^2 / a^2 ...                        
                     * ( ( Wn(i, j + 1) - 2 * Wn(i, j) + Wn(i, j - 1) ) / dt^2 + b * Wn(i - 1, j - 1) ) ...
                     + 2 * Wn(i, j) - Wn(i - 1, j); 
    end

end

figure('Name', 'Numerical solution'); 
surface(t, x, Wn, 'edgecolor', 'none'); 
colormap(jet(256)); 
colorbar; 
xlabel('t'); 
ylabel('x'); 
title('Wn(x, t) - numerical solution');

The plot of the numerical solution is shown here.

The two plotted graphs are not the same, which is proof that I did something wrong in the algorithm. The problem is, I can't find the errors. Please help me find them.

To summarize, please help me change the code so that the two plotted graphs become approximately the same. Thank you for your time.

Answer 1

The finite difference discretization of w_tt = a^2 * w_xx - b*w is

( w(i,j+1) - 2*w(i,j) + w(i,j-1) ) / dt^2
  = a^2 * ( w(i+1,j) - 2*w(i,j) + w(i-1,j) ) / dx^2 - b*w(i,j) 

In your order this gives the recursion equation

w(i,j+1) = dt^2 * ( (a/dx)^2 * ( w(i+1,j) - 2*w(i,j) + w(i-1,j) ) - b*w(i,j) )
            +2*w(i,j) - w(i,j-1)

The stability condition is that at least a*dt/dx < 1. For the present parameters this is not satisfied, they give this ratio as 2.6. Increasing the time discretization to 1000 points is sufficient.

Next up is the boundary conditions. Besides the two leading columns for times 0 and dt one also needs to set the values at the boundaries for x=0 and x=1. Copy also them from the exact solution.

Wn(:,1:2) = Wa(:,1:2); 
Wn(1,:)=Wa(1,:);
Wn(end,:)=Wa(end,:);

Then also correct the definition (and use) of b to that in the source

b      = - (lambda^2 * a^2) + mu^2;  

and the resulting numerical image looks identical to the analytical image in the color plot. The difference plot confirms the closeness