MATLAB求解偏微分方程的方法是数值解法，软件自带有求解亥姆霍兹方程的例程，在命令行中输入edit PDEDEMO2查看代码，在该例程里边界条件是一个带方形孔的单位圆。

http://s15/mw690/003dHy9Ozy6P6oIjJ2m6e&690

 

最终得到的结果是：

 

 

 

http://s16/mw690/003dHy9Ozy6P6oJJtIHaf&690

http://s1/mw690/003dHy9Ozy6P6oL8e2Ya0&690

椭圆型偏微分方程的一般形式为

http://s11/mw690/003dHy9Ozy6P6oNgo542a&690

 

使用assempde()函数来求解这类偏微分方程，调用格式为

http://s5/mw690/003dHy9Ozy6P6oQ67VG24&690

 

求解的边界条件由函数b确定，网格类型由p、e和t确定，c、a、f是椭圆型偏微分方程的参数，对于亥姆霍兹方程，c=1，a=-k^2，f=0。

为了方便查看，我将源代码做了注释，并删掉动态图的部分：

%% 在有方形孔的单位圆中求解亥姆霍兹方程

% 用assempde函数解椭圆型偏微分方程

% function in the Partial Differential Toolbox(TM).

%

% The Helmholtz equation, an elliptic equation that is the time-independent

% form of the wave equation, is

%

% 亥姆霍兹方程的形式为‘-Delta^2u-k^2u = 0’.

%

% Solving this equation allows us to compute the waves reflected by a

% square object illuminated by incident waves that are coming from the

% right.

 

%       Copyright 1994-2014 The MathWorks, Inc.

 

%% Problem Definition

% The following variables will define our problem:

%%

% * |g|: A specification function that is used by |initmesh|. For more

% information, please see the documentation pages for |scatterg| and |pdegeom|.

% * |k|, |c|, |a|, |f|: The coefficients and inhomogeneous term.

g = @scatterg; %定义求解域

k = 60; %入射波有60个波数

c = 1; %椭圆型偏微分方程的三个参数

a = -k^2;

f = 0;

 

%% Boundary Conditions

% 绘制边界和确定边界条件

% condition definition.

figure;

pdegplot(g, 'edgeLabels', 'on'); %绘制pde几何图

axis equal %坐标轴的长度设为相等

title 'Geometry With Edge Labels Displayed'; %设置标题

% Create a pde entity for a PDE with a single dependent variable

numberOfPDE = 1;

pb = pde(numberOfPDE);

% Create a geometry entity

pg = pdeGeometryFromEdges(g);

bOuter = pdeBoundaryConditions(pg.Edges(5:8), 'g', 0, 'q', -60i); %定义边界条件

innerBCFunc = @(thePde, loc, state) -exp(-1i*k*loc.x);

bInner = pdeBoundaryConditions(pg.Edges(1:4), 'u', innerBCFunc);

pb.BoundaryConditions = [bOuter bInner];

 

%%设置网格

% We need a fine mesh to resolve the waves. To achieve this, we

% refine the initial mesh twice.

[p,e,t] = initmesh(g); %将区域分解成三角形网络

[p,e,t] = refinemesh(g,p,e,t); %加密一个三角型网络

[p,e,t] = refinemesh(g,p,e,t); %再次加密，形成的三角型网络更为紧密，但处理的数据更多

figure

pdemesh(p,e,t); %绘制pde网格图

axis equal

%% Solve for Complex Amplitude

% The real part of the vector |u| stores an approximation to a real-valued solution of the

% Helmholtz equation.

u = assempde(pb,p,e,t,c,a,f); %用assempde函数求出来的解

%% Plot FEM Solution

figure

pdeplot(p,e,t,'xydata',real(u),'zdata',real(u),'mesh','off');

colormap(cool) %设置颜色区间，有'jet'(默认),'hsv','hot','cool','spring','gray'等

      

    

     这个方法在确定边界条件上有点麻烦，还是建议用图形化的Partial Differential Toolbox来解偏微分方程。

在命令行中输入pdetool打开编辑器：

 

http://s3/mw690/003dHy9Ozy6P6p2LUki52&690

 

在’pde’-‘pdespecification’里选择Elliptic，椭圆型偏微分方程，并设置好参数使之成为亥姆霍兹方程。

 

 

http://s4/mw690/003dHy9Ozy6P6p47NNVc3&690

绘制边界图形：

放置一个小矩形和一个大圆，并将它们的关系改为相减，通过查看‘boundary’-‘boundary mode’得到边界图

http://s14/mw690/003dHy9Ozy6P6p5fZwx8d&690

 

开始设置边界条件，选中小矩形，在’boundary’-‘specify boundary conditions’里设置，选择狄里克莱边界条件（Ay+By'=C，若B=0，A≠0,则称为第一类边界条件或狄里克莱(Dirichlet)条件；B≠0,A=0，称为第二类边界条件或诺依曼(Neumann)条件；A≠0,B≠0，则称为第三类边界条件或洛平(Robin)条件。）设置h=1，r=10，即在这个边界上解u=10。

http://s12/mw690/003dHy9Ozy6P6p7bnf5fb&690

 

选中大圆，设置边界条件为诺依曼条件，g=0,q=-50i。

http://s4/mw690/003dHy9Ozy6P6p8y0Yrf3&690

 

然后绘制网格并加密：

http://s13/mw690/003dHy9Ozy6P6p9zCvWac&690

 

绘出3D图得到

http://s1/mw690/003dHy9Ozy6P6pbiDQI60&690

 

     这样就绘制出了这个亥姆霍兹方程的解。