Helmholtz 1D Finite Difference Approximation using Algebraic Equation

I am trying to approximate Helmholtz's wave equation using algebraic equation.

I think I have made correct algebraic equation, but it does not work; failed to approximate as figure below.

The boundary condition I used, - ux +iwu = 0 on left edge, ux + iwu = 0 on right edge.

clear
close all

%USER_PAR.m
%% Set Parameters
%%---------------------------------------------------------------
global ax bx w v h k

% define constants
nx = 200; freq = 5;
ax = 0; bx = 1;
w = 2 * pi * freq; v = 1.5; 
k = w / v; h = (bx - ax) / nx;

% define true and source equation
% -Uxx - K^2*u = S
syms x
true(x) = exp(1i*w*(x-1)) + exp(-1i*w*x) - 2;
source(x) = diff(diff(true,x));

true_u = matlabFunction(true);
source = matlabFunction(source);

clear x true

% Main 
%USER_PAR

[A, b] = algebraic_system(source, nx);

u1 = A\b; u1 = u1';
X = linspace(ax, bx, nx+1);
U = true_u(X);

E8 = norm(U(:)-u1(:),inf); E2 = norm(U(:)-u1(:),2)/sqrt(nx/4);
    fprintf('  (nx)=(%3d); (L2,L8)-error = (%.3g , %.3g)\n',nx,E2,E8);

  (nx)=(200); (L2,L8)-error = (11.8 , 9.72)

plot(U); hold on; plot(u1);

%------------------------------------------------------------------

function [A, b] = algebraic_system(source, nx)
global w h ax k

% define A, coefficient.
A = zeros(nx+1, nx+1);

A(1,1) = 2 - h^2*k^2 + 2*h*1i*w; A(1,2) = -2;
A(nx+1, nx+1) = 2 - h^2*k^2 + 2*h*1i*w; A(nx+1, nx) = -2;

for i = 2 : nx
    A(i, i-1) = -1;
    A(i, i) = 2 - h^2*k^2;
    A(i, i+1) = -1;
end

b = zeros(nx+1, 1);
for i = 1: nx+1
    b(i) = source(ax + h*(i-1));
end
b = h^2 * b;
%----------------------------------------------------------
end

[Torsten]

Use bvp4c for real and imaginary part if you have difficulties with the discretization.

Answer 1

[Saarthak Gupta]

Hi,

I understand you are trying to find the numerical solution to the Helmholtz 1D equation using finite difference approximation.

If you are looking for an alternate approach, you can use the ‘bvp4c’ solver which is a finite difference code implementing the three-stage Lobatto IIIa formula to determine a numerical solution to a PDE.

Refer to the following code:

clear
close all

%USER_PAR.m
%% Set Parameters
%%---------------------------------------------------------------
global ax bx w v h k

% define constants
nx = 200; freq = 5;
ax = 0; bx = 1;
w = 2 * pi * freq; v = 1.5; 
k = w / v; h = (bx - ax) / nx;

% define true and source equation
% -Uxx - K^2*u = S
syms x
true(x) = exp(1i*w*(x-1)) + exp(-1i*w*x) - 2;
source(x) = diff(diff(true,x));

true_u = matlabFunction(true);
source = matlabFunction(source);

% clear x true

% Main 
%USER_PAR

X = linspace(ax, bx, nx+1);

solinit = bvpinit(X, @mat4init);
sol = bvp4c(@mat4ode, @mat4bc, solinit);

U = true_u(X);
u1=deval(sol,X);

plot(U);
hold on;
plot(u1);
% hold on;
% plot(u1);
%------------------------------------------------------------------

function dydx = mat4ode(x,y)
global k;
dydx = [y(2)
    -(k.^2)*y(1)];
end

function res = mat4bc(ya,yb)
global w;
res = [-ya(2)+1i*w*ya(1)
    yb(2)+1i*w*yb(1)];
end

function yinit = mat4init(x)
global w;
yinit = [exp(1i*w*(x-1))+exp(-1i*w*x)-2
    -w*exp(-w*x*1i)*1i+w*exp(w*(x-1)*1i)*1i];
end

‘mat4ode’ defines a system of first order PDEs, ‘mat4bc’ defines the boundary conditions, and ‘mat4init’ defines the initial guess for the eigenfunction (u). The solution obtained (‘sol’) can be evaluated on a given interval ([ax,bx], in your case) using ‘deval’.

Please refer to the following MATLAB documentation for more details:

• https://in.mathworks.com/help/matlab/ref/bvp4c.html
• https://in.mathworks.com/help/matlab/boundary-value-problems.html

[Dyuman Joshi]

This answer can be improved -

> Pass the variables to the functions directly, instead of using global

> Not using an inbuilt MATLAB function as a variable - true in this case

%clear
%close all

%USER_PAR.m
%% Set Parameters
%%---------------------------------------------------------------

% define constants
nx = 200; freq = 5;
ax = 0; bx = 1;
w = 2 * pi * freq; v = 1.5; 
k = w / v; h = (bx - ax) / nx;

% define true and source equation
% -Uxx - K^2*u = S
syms x
true0(x) = exp(1i*w*(x-1)) + exp(-1i*w*x) - 2;
%Use the functionality of diff() instead of just using diff() repeatedly
source(x) = diff(true0,x,2);

true_u = matlabFunction(true0);
source = matlabFunction(source);

X = linspace(ax, bx, nx+1);

solinit = bvpinit(X, @(x) mat4init(x,w));
sol = bvp4c(@(x,y) mat4ode(x,y,k), @(ya,yb) mat4bc(ya,yb,w), solinit);

U = true_u(X);
u1=deval(sol,X);

plot(U);
hold on;
plot(u1);
% hold on;
% plot(u1);
%------------------------------------------------------------------

function dydx = mat4ode(x,y,k)
dydx = [y(2)
    -(k.^2)*y(1)];
end

function res = mat4bc(ya,yb,w)
res = [-ya(2)+1i*w*ya(1)
    yb(2)+1i*w*yb(1)];
end

function yinit = mat4init(x,w)
yinit = [exp(1i*w*(x-1))+exp(-1i*w*x)-2
    -w*exp(-w*x*1i)*1i+w*exp(w*(x-1)*1i)*1i];
end