# Helmholtz 1D Finite Difference Approximation using Algebraic Equation

5 views (last 30 days)
Minjae Cho on 28 Aug 2022
Commented: Torsten on 28 Aug 2022
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 on 28 Aug 2022
Use bvp4c for real and imaginary part if you have difficulties with the discretization.