第5章 MATLAB符号工具箱及其应用


MATLAB可以快速进行各种数值运算，尤其是矩阵运算。除此之外，MATLAB还具有符号计算的功能。利用MATLAB符号工具箱，可以对函数进行微分、偏微分和极限计算，定积分和不定积分计算； 还可以求解代数方程和微分方程，对函数做积分变换等。这些功能能够将部分烦琐的推演工作交由计算机进行，从而使人都能够专注于问题的本身，因此具有重要的应用。
5.1MATLAB符号工具箱简介
本节首先分类介绍符号工具箱相关的函数及其应用，然后结合具体的电磁问题进行符号工具的讲解。
5.1.1基本操作命令
本小节给出符号工具箱中最基本的函数并加以介绍。
1. syms创建符号变量和函数
其基本格式如下： 

syms var1 var2 … varN

syms f(var1,var2，…，varN)

例如，下面的命令就创建了两个符号变量x和y。

syms x y

如果要创建一个函数，但不需要给定具体的表达式，则可以用如下形式的定义： 

syms f(x,y)

也可以用如下的方式定义函数，并给定函数的具体表达式： 

syms x y

f(x,y)=x+y


2. sym创建符号变量或者符号数字
例如，可以使用下面的命令创建符号数字2和2/5。

sym('2');sym('2/5')

利用x=sym('x')命令则可以定义符号变量x。





5.1.2表达式化简和替换
本小节给出符号工具箱中的表达式化简和替换用的函数并加以介绍。
1.  pretty将符号表达式以美观的形式呈现出来
基本格式： 

pretty(x)

其中，x为一个符号表达式。

syms a b c

x1=(-b+sqrt(b^2-4*a*c))/(2*a);

x2=(-b-sqrt(b^2-4*a*c))/(2*a);

pretty(x1); pretty(x2);

则显示结果如下： 


-b-(b2-4ac)1/2

--------------------

 2 a



-b+(b2-4ac)1/2

--------------------

 2 a

可见，使用pretty呈现的结果更容易为人所识别。上面的结果就是大家所熟知的一元二次方程的求根公式。
2.  simplify将表达式化简
基本格式： 

simplify(S)

例如： 

syms x

y=(cos(x))^2+(sin(x))^2;

simplify(y)

结果为1。
3.  simple寻找最简表达式
基本格式： 

simple(S)

simple利用内置的各种算法，尝试将表达式化简，并且将最简短的形式呈现出来。因此，是一个强有力的化简工具。
4.  subs替换命令
基本格式： 

subs(s,old,new)

在表达式s中，用new表达式替换全部的old表达式，并重新计算表达式，返回s的一个副本。注意，s本身并不发生变化。

或

subs(s)

将表达式s中的所有符号变量，根据上下文中的定义，或者用工作区中的数值替换，代入到表达式s中，并进行计算。没有赋值的符号变量，依然看作变量。注意，替换之后，s本身并不发生变化。
例如： 

syms x

f=2*x^2-3*x+1;

subs(f,1/3)

ans=2/9

经过上述替换，f的表达式依旧是2*x^2-3*x+1。
还可以替换一个多变量表达式中某个特定变量，如下面例子示意： 

syms x y

f=x^2*y+5*x*sqrt(y);

subs(f,x,3);

ans=9*y+15*y^(1/2)

5. subexpr公共表达式替换命令
基本格式： 

［Y,sigma］=subexpr(X,'sigma')

将表达式X中的公共表达式，用符号变量sigma表示； 简化之后的新的表达式，存放在符号变量Y中。
例如： 

t=solve('a*x^3+b*x^2+c*x+d=0');%求一元三次方程的根

［r,s］=subexpr(t,'s')%对结果进行化简，提取公共表达式，用s表示

结果如下： 

r =

s^(1/3)-b/(3*a)-(- b^2/(9*a^2)+c/(3*a))/s^(1/3)

(- b^2/(9*a^2)+c/(3*a))/(2*s^(1/3))-s^(1/3)/2+(3^(1/2)*(s^(1/3)+

(- b^2/(9*a^2)+c/(3*a))/s^(1/3))*i)/2-b/(3*a)

(- b^2/(9*a^2)+c/(3*a))/(2*s^(1/3))-s^(1/3)/2-(3^(1/2)*(s^(1/3)+

(- b^2/(9*a^2)+c/(3*a))/s^(1/3))*i)/2-b/(3*a)



s =

((d/(2*a)+b^3/(27*a^3)-(b*c)/(6*a^2))^2+(- b^2/(9*a^2)+c/(3*a))^3)^(1/2)-

b^3/(27*a^3)-d/(2*a)+(b*c)/(6*a^2)

作为对比，可以在MATLAB下显示t的表达式，可以发现： 通过将三个根中的公用表达式用s表示出来，大大简化了t的形式。
5.1.3微积分运算
本小节给出符号工具箱中的微积分运算函数并加以介绍。
1. diff微分运算
基本格式： 
diff(expr)，对表达式和默认的变量做导数运算。
diff(expr,n)，对表达式及默认的变量做n阶导数运算。
diff(expr,var,n)，对表达式及指定的变量做n阶导数运算。
例如： 

syms x n

f=besselj(n,x);

diff(f)

ans=

(n*besselj(n,x))/x-besselj(n+1,x)

上面的运算给出了贝塞尔函数的导数形式，这个结果在实际应用中也非常重要！例如，它可以把贝塞尔函数导数的根的问题转化为贝塞尔函数的根的问题，从而简化求解。
下面的例子给出了偏微分运算的情况。

syms x y

f=sin(x)^2+cos(y)^2;

diff(f)

此时，因为没有给定求导变量，MATLAB默认选择距离字母x最近的一个变量进行求导。

ans=

2*cos(x)*sin(x)

也可以指定变量做偏微分运算。例如： 

syms x y

f=sin(x)^2+cos(y)^2;

diff(f,y)

则运行结果如下： 

-2*cos(y)*sin(y)

如果要做高阶导数运算，只需要指定相应的变量，并给出阶数即可。例如： 

syms x y

f=sin(x)^2+cos(y)^2;

diff(f,y,2)

ans=

2*sin(y)^2-2*cos(y)^2

2. int不定积分和定积分运算
基本格式： 
int(expr,var)，对指定的表达式，针对指定的变量做不定积分。
int(expr,var,a,b)，对指定的表达式，针对指定的变量，做积分限为a和b的定积分。
例如： 

syms x

f=sin(x)^2;

int(f)

ans=

x/2-sin(2*x)/4

再如： 

syms x y n

f=x^n+y^n;

int(f)

ans=x*y^n+(x*x^n)/(n+1)

MATLAB根据距离字母x的距离，选择积分变量。也可以指定相应的变量，做不定积分。例如： 

syms x y n

f=x^n+y^n;

int(f,y)

ans=x^n*y+(y*y^n)/(n+1)

如果进行定积分，只需将积分上下限传入到函数int的最后两个参数即可。

syms x y n

f=x^n+y^n;

int(f,1,10);

5.1.4方程求解
本小节给出符号工具箱中用来进行方程求解的函数并加以介绍。
1. solve求解代数方程
基本格式： 
S=solve(eqn,var)，将var当作未知数，求解方程，并将结果赋予变量S。
在只有一个符号变量的情况下，可以不注明未知数var。例如： 

syms x

solve(x^3-6*x^2==6-11*x);

注意： MATLAB使用“==”来定义一个方程，然后即可以用solve进行求解。如果没有给出方程右侧的表达式，solve认为右侧即为0。例如： 

syms x

solve(x^3-6*x^2-6);

ans=

1

2

3

当方程中有多个符号变量时，必须指定对哪个变量进行求解。例如： 

syms x y

solve(6*x^2-6*x^2*y+x*y^2-x*y+y^3-y^2==0,y);

ans=

1

2*x

-3*x


2. solve求解代数方程或方程组
基本格式： 

s=solve(eqns,vars)

syms x y z

［x y z］=solve(z==4*x,x==y,z==x^2+y^2)

如果要改变solve返回的根的顺序，则可以在方程组后面把变量顺序指定。例如： 

syms x y z

［y z x］=solve(z==4*x,x==y,z==x^2+y^2,y,z,x)

3. dsolve常微分方程求解
基本格式： 
dsolve(eq,cond)，对微分方程按照给定的初始（边界）条件进行求解。如果没有给定初始条件，则给出方程的通解形式。
例如： 

syms y(t)

dsolve(diff(y)==y)

ans =

C5*exp(t)

下面的代码给出了初始条件： 

syms x(s) a

x=dsolve(diff(x)==-a*x,x(0)==1)

x=

exp(-a*s)

dsolve还可以求解常微分方程组，如下示意。该代码将求解之后的结果放入一个结构S中，并显示出来。

syms f(t) g(t)

S=dsolve(diff(f)==f+g,diff(g)==-f+g,f(0)==1,g(0)==2);

［S.g;S.f］

输出结果如下： 

ans=

2*exp(t)*cos(t)-exp(t)*sin(t)

exp(t)*cos(t)+2*exp(t)*sin(t)

5.1.5特殊函数
MATLAB符号工具箱内置许多特殊函数，现将最基本的函数或与电磁理论相关的部分函数列举如下。
1.  特殊函数列表mfunlist
用于显示所有的特殊函数列表。可以通过此命令了解MATLAB所包含的所有特殊函数。如下即给出了MATLAB的输出结果。其中，第一列表示函数名称，第二列表示参数，第三列是英文的函数解释。

bernoullin Bernoulli Numbers

bernoullin,zBernoulli Polynomials

BesselIx1,xBessel Function of the First Kind

BesselJx1,xBessel Function of the First Kind

BesselKx1,xBessel Function of the Second Kind

BesselYx1,xBessel Function of the Second Kind

Betaz1,z2Beta Function

binomialx1,x2Binomial Coefficients

EllipticF -z,kIncomplete Elliptic Integral,First Kind

EllipticK -kComplete Elliptic Integral,First Kind

EllipticCK -kComplementary Complete Integral,First Kind

EllipticE -kComplete Elliptic Integrals,Second Kind

EllipticE -z,kIncomplete Elliptic Integrals,Second Kind

EllipticCE -kComplementary Complete Elliptic Integral,Second Kind

EllipticPi -nu,kComplete Elliptic Integrals,Third Kind

EllipticPi -z,nu,kIncomplete Elliptic Integrals,Third Kind

EllipticCPi -nu,kComplementary Complete Elliptic Integral,Third Kind

erfczComplementary Error Function

erfcn,zComplementary Error Functions Iterated Integrals

CizCosine Integral

dawsonxDawsons Integral

PsizDigamma Function

dilogxDilogarithm Integral

erfzError Function

eulernEuler Numbers

eulern,zEuler Polynomials

EixExponential Integral

Ein,zExponential Integral

FresnelCxFresnel Cosine Integral

FresnelSxFresnel Sine Integral

GAMMAzGamma Function

harmonicnHarmonic Function

ChizHyperbolic Cosine Integral

ShizHyperbolic Sine Integral

GAMMAz1,z2Incomplete Gamma Function

Ln,xLaguerre

Ln,x1,xGeneralized Laguerre

WzLamberts W Function

Wn,zLamberts W Function

lnGAMMAzLogarithm of the Gamma function

LixLogarithmic Integral

Psin,zPolygamma Function

SsizShifted Sine Integral

SizSine Integral

Zetaz(Riemann) Zeta Function

Zetan,z(Riemann) Zeta Function

Zetan,z,x(Riemann) Zeta Function



Orthogonal Polynomials

Tn,xChebyshev of the First Kind

Un,xChebyshev of the Second Kind

Gn,x1,xGegenbauer

Hn,xHermite

Pn,x1,x2,xJacobi

Pn,xLegendre

需要指出的是，符号工具箱所提供的函数，其自变量可以是符号变量或者表达式。这个是与普通函数不同的。
2.  特殊函数计算mfun
语法： 

mfun(function,para1,para2,…，paran)

例如： mfun('besselj',0,0)，返回值为1，表明J0(0)=1。
3. 正交多项式
在MATLAB中，内置了多种正交多项式。其格式如下： 
 车贝雪夫第一类、第二类多项式： T(n,x),U(n,x)
 厄米特多项式： H(n,x)
 拉盖尔多项式： L(n,x)
 广义拉盖尔多项式： L(n,a,x)
 勒让德多项式： P(n,x)
在上面的公式中，n表示整数，一般是多项式的阶数； x是一个实数； 广义拉盖尔多项式中的a是一个常数。
5.1.6绘制符号函数的图像
本小节给出符号工具箱中专门用于绘制符号函数图像的绘图函数，并对其应用加以介绍。
1.  ezplot绘制曲线
(1) ezplot(f,［xmin xmax］)，在给定的自变量的区间绘制函数f。如果不指定区间，默认的区间是［-2π，2π］。

syms x

ezplot(x^3-6*x^2+11*x-6);

图51(a)展示了MATLAB绘制的曲线效果。
(2) ezplot(f,［xmin,xmax,ymin,ymax］)，在给定的变量的区间绘制f对应的隐函数。
例如： 

syms x y

ezplot((x^2+y^2)^4==(x^2-y^2)^2,［-1 1］)

图51(b)展示了MATLAB绘制的曲线效果。


图51ezplot绘制的曲线举例


(3) ezplot(x,y,［tmin,tmax］)，在给定的参变量的区间绘制曲线。其中，x、y均为t的函数。
例如： 

syms t

x=t*sin(5*t);

y=t*cos(5*t);

ezplot(x,y,［0 5］);

图52(a)给出了用ezplot具体绘制的曲线形状。
2.  ezplot3绘制三维曲线
基本格式： 
ezplot3(x,y,z,［tmin,tmax］)，绘制参变量t确定的三维空间曲线。
例如： 

syms t

ezplot3(t^2*sin(10*t),t^2*cos(10*t),t );

图52(b)给出了用ezplot3绘制的曲线。


图52含参曲线绘制举例


3.  ezsurf绘制曲面
基本格式： 
(1) ezsurf(f,domain)，在指定的区域内绘制函数f对应的曲面。默认的区域是-2π<x<2π，-2π<y<2π。

syms x y

ezsurf(x^2+y^2)

图53(a)就是ezsurf绘制的一个曲面示意。
(2) ezsurf(x,y,z,［smin,smax,tmin,tmax］)，绘制由参变量s和t确定的曲面。

除了上述函数之外，MATLAB符号工具箱还提供了ezcontour、ezcontourc、ezmesh、ezmeshc、ezploar、ezsurfc等函数，供绘图使用。其基本格式与ezsurf等相似，应用起来非常方便。图53(b)就是ezmeshc绘制的一个曲面。


图53曲面绘制举例


5.2变换电磁理论
如前所述，光学变换实际上是一个广义的坐标变换。它牵涉两个空间，如图54所示。一个是虚拟空间，也称为电磁空间，一般认为是真空(也可以包含介质)，用X=(x，y，z)表示； 另外一个是物理空间，也称为材料空间，用X′=(x′，y′，z′)表示。光学变换的本质就是将虚拟空间的部分区域（如图54（a）中黑色实线包围的区域）映射到物理空间的另外一个区域（如图54（b）中黑色实线与灰色实线之间的部分），由于麦克斯韦方程组的协变性，此空间的映射可以用对应区域的材料电磁参数的改变来表示。从电磁理论的角度来看，物理空间中填充有适当介质的问题，与虚拟空间无介质填充的电磁问题完全等价。这样，人们就可以研究虚拟空间的电磁现象，而通过坐标变换将其变换到物理空间，并通过实际人工电磁材料来实现。在这个过程中，物理本质不变。这就是基于光学变换理论设计新型电磁器件的思路。假设光学变换前空间的电磁参数为
ε=、
μ=
，则变换后物理空间的电磁参数为


ε=′=Aε=
ATdet(A)(51)


μ=的变换公式同上。其中，A=X′X=(x′，y′，z′)(x，y，z)为雅克比矩阵。


图54变换电磁理论牵涉的虚拟空间和物理空间


当虚拟空间与物理空间都采用直角坐标系时，采用式(51)比较方便，直接可以得到物理空间的材料参数。但实际应用中，由于具体边界条件的限制，为方便起见，人们往往使用正交坐标系或者曲线坐标系，材料的参数也都使用相应坐标系下面的表示形式。因此，在使用光学变换的过程中，往往牵涉许多直角坐标系与曲线坐标系下的张量转换问题，比较烦琐。

假设附加在X和X′系上的正交坐标系为R=(u1，u2，u3)与R′=(u′1，u′2，u′3)，它们与直角坐标系单位基矢之间的转换矩阵为


T=Λ(1/hi)XRT(52)


其中： Λ表示对角阵。
因此有


［e1e2e3］T=T［exeyez］T


hi为对应于该正交坐标系的拉梅系数，有


hi=xui2+yui2+zui2i=1,2,3(53)


则由链式求导法则可得


A=X′R′R′RRX(54)

det(A)=detX′R′detR′RdetRX(55)


将上述公式代入式(51)，并考虑到正交坐标系与直角坐标系之间基矢转换矩阵的正交性，则可以得到变换后的正交系下，材料的电磁参数表示形式为


ε=′ort=BRXε=RXTBTdetX′R′detR′RdetXR(56)


若虚拟空间填充的介质为各向同性，则上式还可以进一步简化为


ε=′ort=εBΛ(1/h2i)BTdetX′R′detR′RdetXR(57)


其中


B=Λ(h′i)R′R


于是，材料在直角坐标系下的张量分量就可以通过坐标系的旋转而得到，即


ε=′=T′Tε=′ortT′(58)


利用式(58)，可以得到直角坐标系下的分量表示； 利用式(59)，可以得到正交坐标系下的分量表示。一般曲线系的推导，与此类似。需要注意的是，具体使用时，加撇或者不加撇都适用，即公式对虚拟空间和物理空间都是有效的。


ε=′ort=T′ε=′T′T(59)


5.3基于符号工具箱的变换电磁理论推演

如前所述，MATLAB具有功能强大的符号计算功能，能够进行如微积分、线性代数、积分变换、化简等操作。由于在光学变换中牵涉众多张量转换的操作，大多通过矩阵来进行，因此，利用MATLAB的符号工具箱，将会极大提高推算速度和准确度，对利用复杂变换设计新型电磁器件具有重要的意义。
5.3.1正交坐标系与直角坐标系下材料张量的转换

作为一个简单的例子，对式(58)在柱坐标系的具体情况进行简单的推导。假设已经知道了柱坐标系下的材料表达式为 
ε==Λ(εr，εφ，εz)，考虑用符号推演的方法得到其在直角坐标系下的表达式。对应的MATLAB程序段如下： 

syms x y z rho phi%坐标变量定义

syms er ep ez%材料参数定义

x=rho*cos(phi);%变量关系定义

y=rho*sin(phi);

R=［rho phi z］;

X=［x; y; z］;

T=jacobian(X,R);%雅克比矩阵计算

H=diag(［1 1/rho 1］);%拉梅系数矩阵

T=H*T.';

epsilon=diag(［er ep ez］);

epsilon=T.'*epsilon*T;%推算结果

simple(epsilon)%结果化简

式(510)就是根据MATLAB的输出表达式写出的具体结果。



cos2θεr+sin2θεθcosθsinθ(εr-εθ)0

cosθsinθ(εr-εθ)sin2θεr+cos2θεθ0

00εz(510)



可以看出，这个过程快速、准确。光学变换过程中牵涉众多此类变换，因此，使用符号工具箱，将极大地提高推算效率，使研究者有更多的时间和精力投入到变换本身的优化中去，而不是拘泥于烦琐的数学推导过程。当使用一般曲线坐标系时，这个优势更为突出。

类似的，假设已经知道了球坐标系下的材料表达式为 ε==Λ(εr，εθ，εφ)，也可以用符号推演的方法得到其在直角坐标系下的表达式。对应的MATLAB程序段如下： 

syms x y z r phi theta%坐标变量定义

syms er et ep%材料参数定义

x=r*sin(theta)*cos(phi);%变量关系定义，x

y=r*sin(theta)*sin(phi);%变量关系定义，y

z=r*cos(theta);%变量关系定义，z

R=［r theta phi］;

X=［x; y; z］;

T=jacobian(X,R);%雅克比矩阵计算

H=diag(［1 1/r 1/(r*sin(theta))］);%拉梅系数矩阵

T=H*T.';

epsilon=diag(［er et ep］);

epsilon=T.'*epsilon*T;%推算结果

simple(epsilon)%结果化简

式(511)就是根据MATLAB的输出表达式写出的具体结果。



sin2θcos2φεr+cos2θcos2φεθ+sin2φεφcosφsinφ(sin2θεr+cos2θεθ-εφ)

sinθcosθcosφ(εr-εθ)

cosφsinφ(sin2θεr+cos2θεθ-εφ)sin2θsin2φεr+cos2θsin2φεθ+cos2φεφ

sinθcosθsinφ(εr-εθ)

sinθcosθcosφ(εr-εθ)sinθcosθsinφ(εr-εθ)

cos2θεr+sin2θεθ
(511)



5.3.2变换电磁理论的符号推演
本节采用MATLAB符号工具箱，对几种典型的变换光学理论设计的电磁装置进行电磁参数的推演。为方便起见，仅考虑介电常数。磁导率的考虑与之相似。
1. 柱状隐形装置

首先考虑Pendry等所设计的柱形隐形大衣。采用柱坐标系的表示，隐形大衣的内外半径分别是
a和b，对应的空间映射函数为


ρ′=b-abρ+a，φ′=φ，z′=z



根据式(57)，并在MATLAB下进行编程，得到柱坐标系下的材料电磁参数为


Λ(ρ′-a)ρ′，ρ′ρ′-a，(ρ′-a)b2ρ′(b-a)2ε(512)


这个公式与文献中的结果完全一致。

完成上述符号计算的MATLAB代码如下： 

syms x y z x1 y1 z1 rho phi rho1 phi1 eps a b c%坐标变量定义

x=rho*cos(phi); x1=rho1*cos(phi1);%直角坐标系与柱坐标系间的关系

y=rho*sin(phi); y1=rho1*sin(phi1);%直角坐标系与柱坐标系间的关系

R=［rho phi z］;R1=［rho1 phi1 z1］;%定义几个向量

X=［x; y; z］;X1=［x1; y1; z1］;

det1=det(jacobian(X1,R1));%计算行列式的值

rho1=rho*(b-a)/b+a; phi1=phi; z1=z;%虚拟空间与物理空间的变换关系

R1=［rho1 phi1 z1］;

T=jacobian(X,R);%雅克比矩阵计算

H=diag(［1 rho 1］);%拉梅系数矩阵

H1=diag(［1 rho1 1］);%拉梅系数矩阵

B=H1*jacobian(R1,R);%计算矩阵B

det2=det(jacobian(R1,R));%计算行列式

eps1=eps*B*inv(H*H)*B.'*det(jacobian(X,R))/det1/det2;%光学变换的公式

clear phi1 z1 rho1%清除物理空间的变量

syms phi1 z1 rho1%重新定义这些变量

phi=phi1;z=z1;rho=(rho1-a)*b/(b-a);%虚拟空间与物理空间的变换关系

simple(subs(eps1))%结果化简，并且将最终结果用物理空间的变量表示出来

2. 柱坐标系下的复杂变换
考虑另外一个比较复杂的情况，以验证方法的通用性。在柱坐标系下，采用的变换函数为


ρ′=f(ρ)

φ′=g(φ)

z′=h(z)(513)


下面是MATLAB代码，用于推导物理空间的电磁参数：  

syms x y z x1 y1 z1 rho phi rho1 phi1 eps a b c%坐标变量定义

x=rho*cos(phi); x1=rho1*cos(phi1);%直角坐标系与柱坐标系间的关系

y=rho*sin(phi); y1=rho1*sin(phi1);%直角坐标系与柱坐标系间的关系

R=［rho phi z］;R1=［rho1 phi1 z1］;%定义向量

X=［x; y; z］;X1=［x1; y1; z1］;

det1=det(jacobian(X1,R1));%计算行列式的值

rho1=sym('f(rho)');%虚拟空间与物理空间的变换函数设置

phi1=sym('g(phi)');

z1=sym('h(z)');

R1=［rho1 phi1 z1］;%再次给R1赋值，包含变换函数

T=jacobian(X,R);%雅克比矩阵计算

H=diag(［1 rho 1］);%拉梅系数矩阵

H1=diag(［1 rho1 1］);%拉梅系数矩阵

B=H1*jacobian(R1,R);%定义B矩阵

det2=det(jacobian(R1,R));

eps1=eps*B*inv(H*H)*B.'*det(jacobian(X,R))/det1/det2;%计算电磁参数

simple(eps1)%简化表达式

根据运行结果整理的材料参数为


ε=′=Λ(f′2(fg′)2/ρ2h′2)ερff′g′h′(514)


式(514)与理论推导结果相一致。将其中的不加撇变量用加撇变量表示出来，即可得到物理空间的具体表达式，可以直接指导设计和后续仿真。
3. 直角坐标系下的二维复杂映射
考虑直角坐标系下，二维情况中最一般的映射，也即虚拟空间和物理空间按照如下函数进行变换




x′=f(x，y)

y′=g(x，y)

z′=z(515)



那么根据变换光学理论，容易得到物理空间中的介电常数和磁导率。下面的MATLAB代码就可以完成上述功能。

syms x y z x1 y1 z1 eps%坐标变量定义

x1=sym('f(x,y)');
%定义变换函数

y1=sym('g(x,y)');

z1=z;

A=jacobian(［x1 y1 z1］,［x y z］);%计算雅克比矩阵

eps1=eps*A*A.'/(det(A));%应用式（5-1）

simple(eps1)%化简表达式


根据MATLAB运行结果，写出介电常数张量为



εf2x+f2ygyfx-gxfygxfx+gyfygyfx-gxfy0

gxfx+gyfygyfx-gxfyg2x+g2ygyfx-gxfy0

001gyfx-gxfy(516)


上述结果与理论推导完全一致。

本节内容采用MATLAB下的符号工具箱对光学变换所涉及的各个转换公式进行辅助推演。解析对比和数值仿真验证都证实了上述方法的有效性。该方法快速、高效、准确，可以将科研人员从烦琐的数学推导中解放出来。本节的结果对于应用复杂光学变换、光学变换函数的优化和新型电磁器件的设计有重要的指导作用。
5.3.3介电常数张量的对角化
由变换光学的理论可知，其生成的材料都是非均匀各向异性的磁性材料，需要用人工电磁材料的方法实现。由式(51)可以看出，介电常数张量为对称矩阵，因此，在具体实现的时候，可以利用本征值问题将其在主轴坐标系中对角化。由于只需考虑三个方向的电磁响应特性，这将大大降低材料设计的难度，有利于加工人工电磁材料。
1. 一般情况下的张量对角化
假设光学变换理论中得到的材料的电磁参数为一满阵形式(对称阵)，即


ε==
εxxεxyεxz

εxyεyyεyz

εxzεyzεzz(517)



利用下面的MATLAB代码可以计算该矩阵的本征值，并得到对应的本征向量。由于表达式过长，所以不再列出。

syms epxx epxy epxz epyy epyz epzz%定义符号变量

eps=［epxx epxy epxz; epxy epyy epyz; epxz epyz epzz］;%式（5-17）

［V D］=eig(eps);%计算本征值和本征向量

上面代码中，V中的列向量给出的就是本征向量； D中的对角线元素就是相应的本征值； 以各个本征向量(列向量)为基所构成的新的坐标系就是主轴坐标系。具体在使用的时候，还可以将各个列向量做归一化处理。
2. 直角坐标系中二维情况下的张量对角化
5.3.2节中已经给定了直角坐标系中二维情况下最一般的变换光学的结果，也就是式(516)，也可以简写作式(518)。接下来，采用MATLAB符号工具箱找出其特征值和特征向量。如下： 

syms epxx epxy epyy epzz%定义符号变量

eps=［epxx epxy 0; epxy epyy 0; 0 0 epzz］;%式（5-18）

［V D］=eig(eps);%计算本征值和本征函数



ε==
εxxεxy0

εxyεyy0

00εzz(518)


本征值矩阵D为


Λεzz，εxx+εyy-(εxx-εyy)2+4ε2xy2，εxx+εyy+(εxx-εyy)2+4ε2xy2(519)


对应的本征向量矩阵V为


0εxx-εyy-(εxx-εyy)2+4ε2xy2εxyεxx-εyy+(εxx-εyy)2+4ε2xy2εxy

011

100(520)



如前所述，以式(520)各个列向量为基所构成的坐标系就是主轴坐标系。在这个坐标系下，介电常数张量最简，如式(519)所示。
5.4椭球坐标系的MATLAB辅助分析

椭球坐标系是正交坐标系中较为复杂的一种，在研究包含椭圆边界的电磁问题时非常重要。下面将利用MATLAB符号工具箱来辅助学习椭球坐标系。
5.4.1椭球坐标系中的坐标平面

解决椭球问题时，一般使用椭球坐标系。在这个坐标系中，空间任意一点可以用三个两两正交的曲面的交点表示出来。这三个曲面在直角坐标系中可以统一表示为


x2a2+u+y2b2+u+z2c2+u=1(521)



其中： a，b，c为已知常数，且假定a>b>c。在x，y，z均为确定数值的时候，式(521)是关于u的方程，其有三个不同的实根，分别为ξ，η，ζ，这三个实根也是点(x，y，z)对应的椭球坐标(ξ，η，ζ)。并且这三个根分别有不同的取值范围，即


ξ≥-c2，-c2≥η≥-b2，-b2≥ζ≥-a2


各自对应的曲面方程单独写作



x2a2+ξ+y2b2+ξ+z2c2+ξ=1

x2a2+η+y2b2+η+z2c2+η=1

x2a2+ζ+y2b2+ζ+z2c2+ζ=1
(522)


ξ，η，ζ分别为常数时对应的平面分别为椭球面、单叶双曲面、双叶双曲面，并且均与下面的椭球面共焦，这个椭球面的表达式为



x2a2+y2b2+z2c2=1(523)



它对应的正是ξ=0的情形。
下面的MATLAB代码在直角坐标系中任意选择了一点(6，6，6)，然后计算并绘制通过这一点的三个曲面。通过对这些曲面的观察，就容易理解椭球坐标系中三个坐标平面的形状和特征。需要指出的是，为了绘制式(522)所对应的各个坐标平面，采用了等值面的绘制方式。

a=7;b=5;c=3;%设置椭球系的三个常数

x=6;y=6;z=6;%任意选择一点（6,6,6）

syms u;%定义符号变量u

S=solve(x^2/(a^2+u)+y^2/(b^2+u)+z^2/(c^2+u)==1,u,'Real',true);

%求解u的三个实数解

S=double(S);%将解由符号形式转化为数值

ksi=S(1);eta=S(2);zeta=S(3);%分别赋予ksi、eta、zeta等三个椭球系下的坐标

f1=@(x,y,z) x.^2/(a^2+ksi)+y.^2/(b^2+ksi)+z.^2/(c^2+ksi)-1;%函数表达式

［x,y,z］=meshgrid(-15:.2:15,-15:.2:15,-15:.2:15);%画图范围

v=f1(x,y,z);%定义三元函数f1

h=patch(isosurface(x,y,z,v,0));%绘制v=0的等值面，即式（5-22）对应坐标平面

set(h,'FaceColor','r','EdgeColor','none');%设置坐标面的属性

xlabel('x');ylabel('y');zlabel('z');%标注坐标轴

hold on;%准备继续绘制其他曲面

%以下用类似的方法，分别绘制eta和zeta等于常数的坐标曲面

f2=@(x,y,z) x.^2/(a^2+eta)+y.^2/(b^2+eta)+z.^2/(c^2+eta)-1;%函数表达式

［x,y,z］=meshgrid(-15:.5:15,-15:.5:15,-15:.5:15);%画图范围

v=f2(x,y,z);

h=patch(isosurface(x,y,z,v,0));

isonormals(x,y,z,v,h)

set(h,'FaceColor','g','EdgeColor','none');

f3=@(x,y,z) x.^2/(a^2+zeta)+y.^2/(b^2+zeta)+z.^2/(c^2+zeta)-1;%函数表达式

［x,y,z］=meshgrid(-15:.2:15,-15:.2:15,-15:.2:15);%画图范围

v=f3(x,y,z);

h=patch(isosurface(x,y,z,v,0));

set(h,'FaceColor','b','EdgeColor','none');

alpha(0.6);%曲面半透明

view(3); axis equal; camlight; lighting gouraud;%其他辅助设置

axis vis3d;

上述程序生成的图像如图55所示。为了方便观察，图中将各个坐标面做了单独呈现，并将最终结果放置在图中。


图55椭球坐标系下的坐标平面及其两两相交的情形




图55(续)



5.4.2与直角坐标系的关系
直角坐标系和椭球坐标系的关系为



x=±(ξ+a2)(η+a2)(ζ+a2)(b2-a2)(c2-a2)

y=±(ξ+b2)(η+b2)(ζ+b2)(c2-b2)(a2-b2)

z=±(ξ+c2)(η+c2)(ζ+c2)(a2-c2)(b2-c2)
(524)



由式(524)可以看出，对于特定的ξ，η，ζ，一共有8个点与其对应。从图55可以看出，这8个点对称分布在8个卦限中。下面的MATLAB代码利用符号计算工具，得到式(524)的结果。它实际上是计算求解式(522)所对应的三元方程组所得到的。

syms a b c x y z ksi eta zeta%定义符号变量

%求解三元方程组，得到x、y、z的表达式

S=solve(x^2/(a^2+ksi)+y^2/(b^2+ksi)+z^2/(c^2+ksi)==1,…

 x^2/(a^2+eta)+y^2/(b^2+eta)+z^2/(c^2+eta)==1,…

 x^2/(a^2+zeta)+y^2/(b^2+zeta)+z^2/(c^2+zeta)==1);

%分别显示x、y、z的表达形式

pretty(S.x(1))%用直观的形式，显示x的一个根，下同

pretty(S.y(1))

pretty(S.z(1))

正如前面所述，该方程组一共有8组根，为方便起见，程序仅仅显示了其中的一组。
5.4.3椭球坐标系中的拉梅系数和拉氏运算

拉梅系数是曲线坐标系中的重要元素，它反映了在空间任意一点，当只有一个坐标变量发生变化时(如ξ，η，ζ)，对应的弧长变化的系数。在椭球坐标系中，拉梅系数为


hξ=(ξ-η)(ξ-ζ)2(ξ+a2)(ξ+b2)(ξ+c2)=(ξ-η)(ξ-ζ)2R(ξ)(525a)

hη=(η-ξ)(η-ζ)2(η+a2)(η+b2)(η+c2)=(η-ξ)(η-ζ)2R(η)(525b)

hζ=(ζ-ξ)(ζ-η)2(ζ+a2)(ζ+b2)(ζ+c2)=(ζ-ξ)(ζ-η)2R(ζ)(525c)


其中


R(s)=(s+a2)(s+b2)(s+c2)(526)



下面的MATLAB代码利用符号工具箱推演了上述公式。

syms a b c x y z ksi eta zeta%定义符号变量

S=solve(x^2/(a^2+ksi)+y^2/(b^2+ksi)+z^2/(c^2+ksi)==1,…

 x^2/(a^2+eta)+y^2/(b^2+eta)+z^2/(c^2+eta)==1,…

 x^2/(a^2+zeta)+y^2/(b^2+zeta)+z^2/(c^2+zeta)==1);%求x、y、z的表达式

x=(S.x(1));%x的表达式

y=(S.y(1));%y的表达式

z=(S.z(1));%z的表达式

R_ksi=sqrt((diff(x,ksi))^2+(diff(y,ksi))^2+(diff(z,ksi))^2);%计算拉梅系数

R_eta=sqrt((diff(x,eta))^2+(diff(y,eta))^2+(diff(z,eta))^2);

R_zeta=sqrt((diff(x,zeta))^2+(diff(y,zeta))^2+(diff(z,zeta))^2);

pretty(simple(R_ksi))%直观显示拉梅系数

pretty(simple(R_eta))

pretty(simple(R_zeta))

知道了拉梅系数，则很容易得到拉普拉斯算符在椭球坐标系下的表达式为


2=4R(ξ)(ξ-η)(ξ-ζ)ξR(ξ)ξ+4R(η)(η-ξ)(η-ζ)ηR(η)η

+4R(ζ)(ζ-ξ)(ζ-η)ζR(ζ)ζ(527)



5.5内部匀质化理论及其MATLAB分析
越来越多的研究和实验表明，在基底材料中混合呈周期性排列且结构尺寸远远小于工作波长的其他材料，可以等效为新的、均质材料，这也是设计和加工新型人工电磁材料的重要方法之一。利用均匀媒质来等效周期性亚波长结构的研究方法被称作等效媒质理论（Effective Medium Theory，EMT）。将混合材料做匀质化处理，有两种方法，即外部匀质化和内部匀质化。外部匀质化把一个非均匀材料向外等效为一种均匀材料； 内部匀质化是针对具有复杂结构的混合材料而言，它是将这些具有复杂构造的“原子”结构等效为具有相同形状的均匀的“颗粒”，等效的前提条件是它们对外部的电磁响应相同。图56给出了两种匀质化方法的差异。



图56内部和外部匀质化示意 


5.5.1双层圆柱等效介电常数的分析
对于双层圆柱的等效介电常数的分析，在准静态近似的情况下，首先利用圆柱坐标系内的分离变量法，将双层圆柱和相同外部尺寸的单个圆柱各处电势计算出来。其次基于内部均质化原理，双层圆柱和相同尺寸的单个圆柱对外电场的响应是等效的，因此可以得到双层圆柱的等效介电常数。
1. 单个介质柱



图57均匀电场中的

介质柱


首先来研究一个无限长圆柱介质置于均匀电场中时介质柱内外的电势分布。如图57所示，设在介电常数为εb的无限大均匀介质中存在电场强度E0，垂直于电场方向放置一根半径为b的无限长直介质圆柱体，其介电常数为ε。下面来求解圆柱介质内外的电势。


根据拉普拉斯方程在柱坐标下分离变量的通解形式，可以得到介质柱内外的电势分布有如下形式



1=-E0r+B1rcosφ

2=A1rcos
(528)



考虑到介质表面的边界条件，当r=b时，1=2且εb1r=ε2r，可以得到如下两个方程



-E0b+B1b=A1b

-εbE0-εbB1b2=εA1
(529)


对式(529)进行求解可得


B1=ε-εbε+εbE0b2＝εr-1εr+1E0b2(530)



其中： εr=εεb，表示相对介电常数。
2. 双层介质柱
如图58所示，设在介电常数为εb的无限大均匀介质中存在电场强度E0，垂直于电场方向放置一根内半径为a、外半径为b的无限长双层圆柱介质，其中内层介电常数为ε1，电势为1； 外层介电常数为ε2，电势为2。

图58均匀电场中双层

圆柱的模型


通过上一部分均匀电场中圆柱介质得到的结论，容易得到如下表达式



1=Arcosφ

2=Br+Crcosφ

0=-E0r+Drcosφ
(531)


依旧根据边界上电势和电位移矢量的法向分量连续，可以得到以下的边界条件： 
(1) 当r=a时，有1=2且有ε11r=ε22r
(2) 当r=b时，有2=0且有ε22r=εb0r
将式(531)分别代入上述两个边界条件，可以得到四个方程，此问题有四个未知数，所以可以得到唯一解。四个联立方程如下



Aa=Ba+C/a

ε1A=ε2(B-C/a2)

Bb+C/b=-E0b+D/b

ε2(B-C/b2)=εb(-E0-D/b2)
(532)



因为方程组(532)是由符号构成的，人工进行计算烦琐而且容易出现错误。利用MATLAB求解符号方程的途径是理想的方法。因此将方程(532)写成矩阵方程形式非常有助于利用MATLAB求解出未知数。将方程组(532)写成矩阵方程的形式为



a-a-1a0

ε1-ε2ε2a20

0b1b-1b

0ε2-ε2b2εbb2A

B

C

D=0

0

-E0b

-E0εb  (533)



内部均质化原理的本质在于，准静态近似下无限长圆柱介质和双层圆柱之间存在一个等量关系，即B1≡D。此时，从外部来看，双层圆柱与单层圆柱的电磁响应是一模一样的，或者说二者是等效的。D的表达式为



D=E0b2b2(ε1+ε2)(ε2-εb)+a2(ε1-ε2)(ε2+εb)b2(ε1+ε2)(ε2+εb)+a2(ε1-ε2)(ε2-εb)(534)


根据式(530)，由于B1≡D，所以很容易得到


εr-１εr+１E0b2=D (535)


为了表示方便和编程简单，设=D/E0b2，则有


εr=1+1-(536)



使用MATLAB编写程序，将所求得的D的值代入式(536)，可以得到最终相对介电常数的值为


εr=εεb=ε2［b2(ε1+ε2)+a2(ε1-ε2)］εb［b2(ε1+ε2)-a2(ε1-ε2)］ (537)


约去式(537)两边重复的εb，ε即为待求解的介电常数，即双层柱的等效介电常数


ε=ε2［b2(ε1+ε2)+a2(ε1-ε2)］b2(ε1+ε2)-a2(ε1-ε2) (538)


如果设内外层圆柱的半径比为rcs=a/b，则式(538)又可以进一步简化为


ε=ε2*［(ε1+ε2)+r2cs(ε1-ε2)］(ε1+ε2)-r2cs(ε1-ε2) (539)



下面的MATLAB代码就实现了上面的代数方程求解并进行化简的过程。

syms a ep_1 ep_2 ep_b b E0 b ep_e

A=［a-a-1/a0;

ep_1 -ep_2ep_2/a^20;

0b1/b-1/b;

0ep_2-ep_2/b^2ep_b/b^2］;%方程系数矩阵

m=［0;0;-E0*b;-ep_b*E0］;%方程右侧部分

C=A＼m;%左除的位置向量

D=C(4);%D的结果

D_m=D/(E0*b^2);%将D化简

［N1,D1］=numden(D_m)%D_m的分子和分母，用于验证

ep_e=(1+D_m)/(1-D_m);%等效介电常数

［N2,D2］=numden(ep_e)%ep_e的分子和分母

3. 双层介质柱等效介电常数的变化曲线
在本部分使用银和二氧化硅两种物质来研究当内外层材料不同时，等效介电常数与内外层圆柱的半径比
rcs的关系。假设波长为405nm，此时银的介电常数为
εm=-4.70+i0.22，二氧化硅的介电常数为εd=2.42。

首先，当外层是银介质，内层是二氧化硅时，等效介电常数可以写作εeff=ε′eff+iε
″eff的形式，在MATLAB中绘制这种情况下等效介电常数的实部ε′eff与虚部ε″eff随着内外层半径比的关系。变化关系如图59所示。

从图59可以看出，这种情况下双层圆柱等效介电常数的虚部接近0，而介电常数实部从-4.70到2.42变化，也就是说外层为金属介质时，等效介电常数的实部不可能超出其组成材料的介电常数的范围。

当内层是银介质，外层是二氧化硅时，等效介电常数也可以写作εeff=ε′eff+iε″eff的形式。同样地，在MATLAB中绘制这种情况下等效介电常数的实部ε′eff与虚部ε″eff随着内外层半径比a/b的关系。变化关系如图510所示。


图59外层为银，内层为二氧化硅时介电

常数的变化趋势




图510内层为银，外层为二氧化硅时介电

常数的变化趋势


从图510可以看出，内层为银介质的情况下，介电常数的变化范围很广，其变化范围不单纯的是在两种材料的介电常数范围内，实际应用时可以根据设计的实际需要来确定内外层材料和内外圆柱半径比。

由上述计算可以看出，采用两种材料就可以搭配出诸多介电常数的组合，从而实现具有特殊功能的电磁器件。因此，内部匀质化方法是实现人工电磁材料的重要方法。
5.5.2双层球结构
与上述做法相似，可以考虑双层介质球的等效介电常数。
1. 单层球介质的电势分布
首先研究在无限大均匀电场中放置的单层介质球的电势分布。如图511所示，设均匀电场的电场强度为
E0，方向沿z轴方向。背景材料介电常数为εb，电势为1； 

图511均匀电场中的

单层介质球
置于均匀电场中的单层介质球的介电常数为ε，电势为2。


利用球坐标系下拉普拉斯方程的通解，各个区域的电势分布为



1=(-E0r+B1r-2)P1(cosθ)


2=A1rP1(cosθ)
(540)


根据当r=b时，1=2且εb1r=ε2r这个条件可以得到如下方程组




-E0b+B1b2=A1b

-εbE0-2εbB1b3=εA1
(541)


求解方程组(541)可以得到


B1=ε-εbε+2εbE0b3＝εr-１εr+2E0b3(542)


其中，εr=εεb，表示相对介电常数。
2. 双层球介质的电势分布
上一部分已经得到了单层球介质的电势分布的表达式，与求解双层柱的方法类似，本部分待求解的是双层球介质每层的电势分布，

图512均匀电场中的

双层介质球
根据单层和双层之间存在的某些等量关系，就可以求得双层球介质的等效介电常数。

如图512所示，设在介电常数为εb的无限大均匀介质中存在电场强度E0，其中放置有一个内半径为a、外半径为b的双层介质球，内层介电常数为ε1，电势为1； 外层介电常数为ε2，电势为2。

根据上一部分得到的结论，可以很容易得到每层的电势表达式为



1=ArP1(cosθ)

2=［Br+Cr-2］P1(cosθ)

0=［-E0r+Dr-2］P1(cosθ)
(543)


根据边界上电势和电位移矢量的法向分量连续，可以得到如下边界条件： 
(1) r=a时，1=2且ε11r=ε22r
(2) r=b时，2=0且ε22r=εb0r
将上述边界条件应用到式(543)中，可以得到如下方程组



Aa=Ba+C/a2

ε1A=ε2(B-2C/a3)

Bb+C/b2=-E0b+D/b2

ε2(B-2C/b3)=εb(-E0-2D/b3)
(544)



为了方便在MATLAB中进行计算，将方程(544)写成矩阵方程，其形式为


a-a-1a20

ε1-ε22ε2a30

0b1b2-1b2

0ε2-2ε2b32εbb3A

B

C

D=0

0

-E0b

-E0εb(545)


利用内部均质化原理可以得知，当双层介质球与同等大小的单层介质球放入同样的均匀电场中时，在准静态近似下引起的响应必然是相同的。所以必有B1≡D。因此可以得到


εr-１εr+2E0b3≡D(546)


为了表示方便，此处设=D/E0b3，可以得到


εr=1+21-=1+2D/E0b31-D/E0b3(547)


利用MATLAB编程求解方程(545)，解出D的表达式为


D=E0b3b3(ε1+2ε2)(ε2-εb)+a3(ε1-ε2)(εb+2ε2)b3(ε2+2εb)(ε1+2ε2)-2a3(ε1-ε2)(ε2-εb)


将D代入式(547)中，有


εr=εεb=ε2［b3(ε1+2ε2)+2a3(ε1-ε2)］εb［b3(ε1+2ε2)-a3(ε1-ε2)］(548)


观察式(548)，约去公式左右两边的εb，可以得到双层椭球的等效介电常数ε为


ε=ε2*［b3(ε1+2ε2)+2a3(ε1-ε2)］［b3(ε1+2ε2)-a3(ε1-ε2)］(549)



下面的MATLAB代码可以实现上述方程组求解和等效介电常数计算的过程。

syms a ep_1 ep_2 ep_b b E0 b ep_e

A=［a-a-1/(a^2)0;

ep_1-ep_22*ep_2/a^30;

0b1/(b^2)-1/(b^2);

0ep_2-2*ep_2/b^32*ep_b/b^3］;%定义系数矩阵

m=［0;0;-E0*b;-ep_b*E0］;%方程右侧部分

C=A＼m;%求未知数向量［A B C D］^T

D=C(4);%D的表达式

D_m=D/(E0*b^3);%简化D

［N1,D1］=numden(D_m)%求D_m的分子和分母

ep_e=(1+2*D_m)/(1-D_m);%式（5-47）

［N2,D2］=numden(ep_e)%求ep_e的分子和分母，以便于验证

5.5.3双层椭球结构的等效介电常数
1. 椭球形介质周围的势函数


接下来，采用类似的方法，考虑共焦双层椭球结构的情形。将双层椭球之间的边界定义为ξ1和ξ2，且注意双层椭球共焦的条件为椭球的三个半轴ai，bi，ci必须满足如下条件


a21-a22=b21-b22=c21-c22=-ξ2(550)


从式(550)可以明显看出ξ1的值为0。
如图513所示，双层椭球和同等大小的单层均质椭球均被置于均匀外电场中。外电场E沿x轴方向。


图513双层椭球及等效的单层均质椭球示意


对于单层椭球而言，设无穷远处电势为∞，很容易得到其表达式为


∞=-Ex=-E(ξ+a2)(η+a2)(ζ+a2)(b2-a2)(c2-a2)(551)


而椭球附近电势可以表示为=∞F(ξ)。朗道所编写的《连续介质的电动力学》一书中已经给出了F(ξ)的最终表达式。F(ξ)为常数或为



F(ξ)=∫∞ξds(s+a2)R(s)(552)


其中


R(s)=(s+a2)(s+b2)(s+c2)


由于受到式(550)所示的椭球共焦条件的限制，所以式(552)中的a，b，c可以替换成任意的ai，bi，ci。此处用a1，b1，c1代替，因此电势的表达式可以写作


=∞C-D2
∫∞ξds(s+a21)R1(s)(553)

R1(s)=(s+a21)(s+b21)(s+c21)(554)



与双层球和双层柱相同，边界上的电势应是连续的，即C-D2∫∞ξds(s+a21)R1(s)在边界上连续。另外，边界上的法向电位移也应该是连续的，电位移矢量D的法向部分表达式为


n·D=-ε1hξξ(555)


其中： hξ为椭球坐标系的度量系数，其表达式为


hξ=(ξ-η)(ξ-ζ)2R1(ξ)(556)


将式(556)代入式(555)可以得到： 要使法向电位移连续，也就是必须使表达式
εC+DR1(ξ)-D2∫∞ξds(s+a21)R1(s)在边界上连续。
2. 单层椭球的电势分布
如图513的右侧所示的单层均质椭球，分析其每层的电势分布。注意到内层椭球F(ξ)只能是常数，否则积分那一部分在原点处存在奇异值。所以根据上述分析可以很容易得到如下电势分布



e=-CEx(ξ≤ξ1)

0=-Ex1-D2∫∞ξds(s+a21)R1(s)(ξ>ξ1)
(557)


根据电势连续和法向电位移连续可以得到以下方程


C=1-D2∫∞ξ1ds(s+a21)R1(s)(558a)

εeC=ε01+DR1(ξ1)-D2
∫∞ξ1ds(s+a21)R1(s)(558b)


上述方程依旧是符号组成的方程，利用MATLAB解方程可以得到D的表达式为


D=2R1(ξ1)(εe-ε0)2ε0+N1R1(ξ1)(εe-ε0)(559)


其中


N1=∫∞ξ1ds(s+a21)R1(s)=2a1b1c1Nx1(560)


其中： Nx1为x方向的退极化因子。
上述求解系数D的过程，可以用MATLAB代码计算如下： 

syms ep_e ep_0 N1R1%定义符号变量

A=［1 N1/2;ep_e ep_0*N1/2-ep_0/R1］;%定义系数矩阵

m=［1;ep_0］;%定义方程右侧的列向量

C=A＼m;%求解代数方程

D=C(2);%D的表达式

pretty(D)%用直观的形式显示系数D

3. 双层椭球的电势分布
图513中左侧双层椭球的电势分布与单层均质椭球的电势分布是类似的，每层的电势为



2=-C2Ex(ξ≤ξ2)

1=-ExC1-D12∫∞ξds(s+a21)R1(s)(ξ2<ξ<ξ1)

0=-Ex1-D02∫∞ξds(s+a21)R1(s)(ξ≥ξ1)
(561)


边界条件仍旧是边界上电势和法向电位移连续，由此可得以下方程组


C2=C1-D12∫∞ξ2ds(s+a21)R1(s)(562a)

ε2C2=ε1C1+D1R1(ξ2)-D12∫∞ξ2ds(s+a21)R1(s)(562b)

C1-D12∫∞ξ1ds(s+a21)R1(s)=1-D02∫∞ξ1ds(s+a21)R1(s)(562c)

ε1C1+D1R1(ξ1)-D12∫∞ξ1ds(s+a21)R1(s)=ε01+D0R1(ξ1)-D02∫∞ξ1ds(s+a21)R1(s)(562d)



观察式(557)与式(561)可以发现，常数D与D0分别代表将均质椭球与双层椭球放置于外加均匀电场中时产生的扰动。当椭球尺寸远小于波长时，在准静态近似情况下应有D≡D0，即二者对椭球外电场的响应相同。这就是内部均质化原理的本质。
利用MATLAB编写程序，解方程组(562)，可以得到D0的表达式，其形式过于复杂，此处暂不给出。再由D≡D0，解得等效介电常数的表达式为



εe=ε11+2R1(ξ2)(ε2-ε1)R1(ξ1)［2ε1+(N1-N2)(ε1-ε2)R1(ξ2)］(563)


其中: N2的值为


N2=∫∞ξ2ds(s+a21)R1(s)=2a2b2c2Nx2


上面已经提到a21-a22=b21-b22=c21-c22=-ξ2且ξ1=0，因此有


R1(ξ1)=a1b1c1，R1(ξ2)=a2b2c2


对式(563)进行化简可得


εe=ε11+(ε2-ε1)λε1+(ε1-ε2)(λNx1-Nx2)(564)


其中


λ=a2b2c2a1b1c1


这就是双层椭球的等效介电常数的最终结果。

下面的MATLAB程序正是对上述过程进行符号推演的具体实现。

syms ep_1 ep_2 ep_0 N1 N2 R1 R2 ep_r ep_e%定义符号变量

A=［-1 1 -N2/2 0;

-ep_2 ep_1 ep_1/R2-ep_1*N2/2 0;

0 1 -N1/2 N1/2;

0 ep_1 ep_1/R1-ep_1*N1/2 ep_0*N1/2-ep_0/R1］;

m=［0;0;1;ep_0］;

C=A＼m;%C=［C2 C1 D1 D0］

D0=C(4);%D0的表达式

A=［1 N1/2;ep_e ep_0*N1/2-ep_0/R1］;

m=［1;ep_0］;

C=A＼m;

D=C(2);%D的表达式

ep_r=solve(D==D0,ep_e);%求解方程，得到ep_e

pretty(simple(ep_r))%显示ep_e的表达式

5.6无限大半平面的衍射——特殊函数的应用
借助于惠更斯菲涅尔原理可以解释和描述光束通过各种形状的障碍物时所产生的衍射现象。
讨论时，通常可以根据光源和考察点到障碍物的距离把衍射现象分为两类。第一类是障碍物到光源和考察点的距离都是有限的，

图514半无限大平板的

衍射示意图
或其中之一为有限的，称为菲涅尔衍射，又称近场衍射； 第二类是障碍物到光源和考察点的距离可认为是无限远的，即实际上使用的是平行光束，这种衍射称为夫琅禾费衍射，又称远场衍射。


采用标量衍射理论，考察平行光垂直入射到一个半无限大平板时发生的衍射情况，如图514所示。平板占据的位置为x′<0，与观察屏之间的距离为d。

则观察屏上电磁波的复振幅分布为


U(x，y)=A∫∞-∞∫∞0e-jkRRdx′dy′(565)


其中： A为半无限大平板上的电磁场的幅度。且有


R=［(x′-x)2+(y′-y)2+d2］1/2


考虑到观察屏到半无限大平面的距离非常大，所以可以将上式近似用牛顿公式表示为


R≈d+(x′-x)2+(y′-y)22d


代入式(565)，并注意到幅度和相位对精度要求的差异，整理得到


U(x，y)=Ade-jkd∫∞-∞∫∞0exp-jπ22λd［(x′-x)2+(y′-y)2］dx′dy′


做变量替换


u=2λd(x′-x)，v=2λd(y′-y)


则有


U(x，y)=Aλ2e-jkd∫∞-ζe-jπ2u2du∫∞-∞e-jπ2v2dv(566)


其中


ζ=2λdx


利用菲涅尔积分的定义，式(566)可以写作


U(x，y)=Aλ2e-jkd∫∞-ζcosπ2u2-jsinπ2u2du∫∞-∞cosπ2v2-jsinπ2v2dv


其中


C(x)=∫x0cosπt22dt，S(x)=∫x0sinπt22dt(567)



考虑到被积函数的奇偶性，即


U(x，y)=Aλ2e-jkdC(∞)-jS(∞)+C(ζ)-jS(ζ)

2C(∞)-j2S(∞)(568)


由于C(∞)=S(∞)=0.5，式(568)可以化简为


U(x，y)=Aλ2e-jkd12-j12+C(ζ)-jS(ζ)(1-j)


所以，观察屏上的光强为


I=A2λ22
12+C(ζ)2+12+S
(ζ)2=I012+C(ζ)2+12+S(ζ)2(569)



下面的MATLAB代码就计算了光场分布的情况。

zeta=-5:0.1:5;%定义zeta区间，它与x成正比

C=mfun('FresnelC',zeta);%计算菲涅尔余弦积分值

S=mfun('FresnelS',zeta);%计算菲涅尔正弦积分值

I0=1;

I=(C+1/2).^2+(S+1/2).^2;%计算光场强度

plot(zeta,I);%绘制曲线

%figure; plot(C,S);%绘制考纽螺线


图515给出了观察屏上x轴上的光场分布。从图515中可以看出，在几何光学中光线照射不到的位置(x<0)，光强不为零； 在靠近几何光学生成的阴影边界处(x>0)，光强是明暗分布的，直到距离边界很远地方，光强才稳定下来。



图515观察屏上x轴方向上看到的情况



事实上，利用符号工具箱中的mfun函数很容易得到
C(∞)=S(∞)=0.5。在MATLAB命令行下输入mfun('FresnelC',inf)或mfun('FresnelS',inf)，则可得到其积分值，即0.5。
另外，当自变量取值一致时，分别以菲涅尔余弦积分和正弦积分作为横纵坐标，可以绘制所谓的考纽螺线，如图516所示。只需在上述代码部分增加一条指令即可(上面代码中注释掉的最后一行)。


图516考纽螺线示意图



5.7小结
本章介绍了MATLAB环境中的符号工具箱。利用该符号工具箱的强大功能，可以辅助进行符号推演，从而将人们从烦琐、易错的劳动中解脱出来，大大提高工作效率。以变换电磁理论的推演、线性方程组求解、张量对角化、椭球坐标系和半无限大平板的衍射为例，详细介绍了各个相关命令、函数及其实际应用。符号推演得到的结果相对复杂，因此要通过各种方式进行简化，并需要人工介入，在应用的时候必须留意。