Helmholtz 1D Finite Difference Approximation using Algebraic Equation
9 次查看(过去 30 天)
显示 更早的评论
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 个评论
Torsten
2022-8-28
Use bvp4c for real and imaginary part if you have difficulties with the discretization.
采纳的回答
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
2023-9-17
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!