Helmholtz 1D Finite Difference Approximation using Algebraic Equation

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

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by