Helmholtz 1D Finite Difference Approximation using Algebraic Equation

7 次查看（过去 30 天）
Minjae Cho
Minjae Cho 2022-8-28
评论： Dyuman Joshi 2023-9-17
Ran in:
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

请先登录，再回答此问题。

采纳的回答

Saarthak Gupta
Saarthak Gupta 2023-9-14
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:
  1 个评论
Dyuman Joshi
Dyuman Joshi 2023-9-17
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by