2

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.

Nikola Ristic
  • 133
  • 4
  • You have a problem with the order of indices. In the declaration of Wn, it is (space,time). In the initialization of the first two time steps, it appears as (time,space). In the update you seem to use (space,time) by the loop headers, but you increment the space index instead of the time index in the update formula. – Lutz Lehmann Oct 02 '22 at 10:52
  • Hello @LutzLehmann. I have changed the double for loop to be: `for i = 2 : (number_of_discrete_space_steps - 1)` `for j = 2 : (number_of_discrete_time_steps - 1)` `Wn(i, j + 1) = dt^2 * ( (a^2 / dx^2) * ( Wn(i + 1, j) - 2 * Wn(i, j) + Wn(i - 1, j) ) - b * Wn(i - 1, j - 1) ) + 2 * Wn(i, j) - Wn(i, j - 1);` `end` `end` but the output graph is still not correct. Do you have any other suggestions? – Nikola Ristic Oct 02 '22 at 12:29
  • The initialization also needs to be changed, `Wn(:,1)=Wa(:,1)`, provided `Wa` is actually in the format (space,time). Are you sure about the term `- b * Wn(i - 1, j - 1)`, it seems to be unsymmetric? /// With slice operations you could eliminate the space loop, this might be faster, especially if you can collect the coefficients so that the operation becomes a convolution. – Lutz Lehmann Oct 02 '22 at 13:08
  • Even after changing the initialization to `Wn(:,1)=Wa(:,1);` `Wn(:,2)=Wa(:,2);` , I still didn't get the correct graph. Also, I still don't know what is a slice operation.I will google it. If you know a good link, please let me know. – Nikola Ristic Oct 02 '22 at 13:24

1 Answers1

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

enter image description here

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51