Helmholtz 1D Finite Difference Approximation using Algebraic Equation
9 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Minjae Cho
el 28 de Ag. de 2022
Comentada: Dyuman Joshi
el 17 de Sept. de 2023
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1109875/image.png)
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);
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 comentario
Torsten
el 28 de Ag. de 2022
Use bvp4c for real and imaginary part if you have difficulties with the discretization.
Respuesta aceptada
Saarthak Gupta
el 14 de Sept. de 2023
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 comentario
Dyuman Joshi
el 17 de Sept. de 2023
This answer can be improved -
> Pass the variables to the functions directly, instead of using global
%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
Más respuestas (0)
Ver también
Categorías
Más información sobre Equation Solving en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!