复变函数与场论简明教程：复变函数与场论的MATLAB求解

复变函数与场论的MATLAB求解

7.1复变函数的MATLAB求解7.2场的MATLAB求解

7.1复变函数的MATLAB求解

1.复数的表示

［例1］将变量a赋值为3+4i，将变量b赋值为10a+5ai。解在命令窗口输入

a=3+4*i（或输入a=3+4i）结果为

a=

3.0000+4.0000i

继续输入

b=10*a+5*a*i

结果为

b=

10.0000+55.0000i

也可以对变量采用模和辐角的方式赋值，即z=rexp(i*θ)，其中r和θ为实数,exp为指数函数。默认变量pi表示圆周率π。［例2］将变量a赋值为10ei(π／4)。

解输入

a=10*exp(i*pi/4)

结果为

a=

7.0711+7.0711i

2.复数的运算

1)实部、虚部、共轭、模和辐角

设z为一复数，MATLAB函数real(z)可以计算z的实部；imag(z)可以计算z的虚部，conj(z)可以计算z的共轭，abs(z)可以计算z的模，angle(z)可以计算z的辐角主值。［例3］求的实部、虚部、共轭、模和辐角主值。

解输入

z=-sqrt(12)-2i；

result=［real(z)；imag(z)；conj(z)；abs(z)；angle(z)］结果为

result=

－3.4641

－2.0000

－3.4641+2.0000i

4.0000

－2.6180

2)加、减、乘、除、幂与方根

设a和b为两复数，n为正整数,则在MATLAB中，a+b可以实现复数加法，a－b可以实现复数减法，a*b可以实现复数乘法，a/b可以实现复数除法,a^n可以计算复数

a的n次幂,a^(1/n)可以计算复数a的n次方根(需要注意的是，该结果仅对应于

在k=0时所对应的根)。［例4］设a=1－2i，b=2+i，求。

解输入

a=1－2i；

b=2+i；

result=(a^3+b^5)/(a^(1/4))

结果为

result=

－48.1545+22.8756i

3)复数方程求根

为了求复数所有的n次方根，可以利用MATLAB函数solve(EQ)。该函数可以对以EQ为表达式的复数方程求解。［例5］求方程z3=1的解。

解输入

EQ=′z^3=1′；

result=solve(EQ)结果为

result=

1

－1/2+1/2*sqrt(－1)*3^(1/2)

－1/2－1/2*sqrt(－1)*3^(1/2)

所得的结果为符号变量，对结果中的符号变量利用函数eval进行数值化，可得三个根。继续输入

nresult=eval(result)

结果为

nresult=

1

－0.5000+0.8660i

－0.5000－0.8660i

该方法相当于求出的三个根。

4)指数、对数和乘幂

MATLAB函数exp(z)可以计算z的指数，log(z)可以计算z以e为底的对数主值，a^b可以计算a的b次幂的主值。

［例6］设a=1－2i，b=5+i，求ea，ln(a+ib)的主值以及ii的主值。

解输入

a=1－2i；

b=5+i；

result=［exp(a)；log(a+i*b)；i^i］结果为

result=

－1.1312－2.4717i

1.0986+1.5708i

0.2079

5)三角函数与反三角函数运算

［例7］设z=1－2i，求sin(z)，sin(z)/cos(z)，tan(z+2π),arctan(tan(z+10π))以及csc(z)的值。

解输入

z=1－2i；

result=［sin(z)；sin(z)/cos(z)；tan(z+2*pi)；

atan(tan(z+10*pi))；csc(z)］结果为

result=

3.1658－1.9596i

0.0338－1.0148i

0.0338－1.0148i

1.0000－2.0000i

0.2284+0.1414i

6)双曲函数与反双曲函数运算

［例8］设z=1－2i，求z的双曲正弦值及

的值。

解输入

z=1-2i；

result=［sinh(z)；(exp(z)-exp(-z))/2］结果为

result=

－0.4891－1.4031i

－0.4891－1.4031i

3.求导运算

利用MATLAB函数syms可以定义符号变量，从而进行求导、积分、求泰勒级数等运算。利用MATLAB函数diff(FUN,n)可以对函数FUN求n阶导数，其中n为正整数。利用MATLAB函数pretty可以将符号类型变量表达为易读的形式。利用MATLAB函数latex可以将符号类型变量表达为LATEX语言形式。［例9］设f(z)=1/z，求f(z)的一阶、二阶、三阶导数。解输入

symsz；

FUN=(1/z)；

result=［diff(FUN,1)；diff(FUN,2)；diff(FUN,3)］结果为

result=

-1/z^2

2/z^3

-6/z^4

继续输入

pretty(diff(FUN))

［例10］设，求f′(z)。

解输入

symsz；

FUN=log(z+sqrt(z^2+1))；

result1=diff(FUN)

结果为

result1=

(1+1/(z^2+1)^(1/2)*z)/(z+(z^2+1)^(1/2))结果为

继续输入

latex(diff(FUN))

结果为

＼left(1+{＼frac{z}{＼sqrt{{z}^{2}+1}}}＼right)＼left(z+＼sqrt{{z}^{2}+1}＼right)^{-1}

即

4.积分运算

利用MATLAB函数int(FUN,z)可以对包含以z为变量的函数FUN求不定积分；函数int(FUN,z,a,b)可以对包含以z为变量的函数FUN求上限为a、下限为b定积分，其中a和b可以为符号变量，也可以为数值变量或数值。［例11］设

解输入

symszab；

FUN=2*z/(1+z^2)^2；

pretty(int(FUN,z))求结果为继续输入

pretty(int(FUN,z,a,b))

结果为继续输入

pretty(int(FUN,z,1,2+i))

结果为

5.泰勒级数展开

［例12］设f(z)=ez，求泰勒级数展开。

解输入

symsz；

FUN=exp(z)；

result=taylor(FUN,6,0)结果为

result=

1+z+1/2*z^2+1/6*z^3+1/24*z^4+1/120*z^5

继续输入

pretty(result)

结果为

6.零点计算

［例13］求p(z)=z3－2z2+3z－4的零点。

解输入

p=［1,－2,3,－4］；

z=roots(p)结果为

z=

1.6506

0.1747+1.5469i

0.1747-1.5469i

7.留数计算

［例14］求

的留数。

解输入

［r,p］=residue(［3,4,3］,［1,-1,1,-1］)结果为

r=

-1.0000-1.0000i

-1.0000+1.0000i

5.0000

p=

0.0000+1.0000i

0.0000-1.0000i

1.0000

即

Res［f(z),i］=－1－i，Res［f(z),－i］=－1+i，Res［f(z),1］=5。［例15］计算积分

其中,C为正向圆周,|z|=2。

解首先计算函数的极点及对应留数。输入

［r,p］=residue(［1,0］,［1,0,0,0,-1］)

结果为

r=

0.2500

0.2500

－0.2500－0.0000i

－0.2500+0.0000i

p=

－1.0000

1.0000

－0.0000+1.0000i

－0.0000－1.0000i

即有四个极点，分别为1，－1，i，－i，它们都在圆周|z|=2内。其对应留数分别为0.25，0.25，－0.25,－0.25。［例16］求的留数。

解输入

symsz；

FUN=1/(z+i)；

r=limit(FUN*(z+i),z,-i)

结果为

r=

1

8.函数图形

［例17］绘制f(z)=z3的图形。

解输入

colormap(hsv(64))

z=cplxgrid(30)；

cplxmap(z,z.^3)

colorbar(′vert′)

输出图形如图7.1所示，其中，f(z)的实部用幅度表示，虚部用颜色表示。图7.1

7.2场的MATLAB求解

1.梯度计算

利用MATLAB函数jacobian(u,x,y,z)可以对以符号化表示的数量场u求其梯度场。若对所得的梯度进行数值化，可得该数量场在确定点处的梯度。［例1］求u(x,y,z)=xy2+yz3的梯度场，并求在点

M(2,－1,1)处的梯度及其沿矢量l=2i+2j－k方向的方向导数。解输入

symsxyz

u=x*y^2+y*z^3；

gradu=jacobian(u,［x,y,z］)结果为

gradu=

［y^2,2*x*y+z^3,3*y*z^2］

即

grad

u=y2i+(2xy+z3)j－3yz2k继续输入

x=2；

y=-1；

z=1；

graduM=eval(gradu)

结果为

graduM=

1-3-3即

gradu|M=i－3j－3k

继续输入

d=graduM*［2,2,－1］′/sqrt(2^2+2^2+(－1)^2)

结果为

d=

－0.3333

即在矢量l=2i+2j－k方向的方向导数为－0.3333。［例2］绘制u(x,y,z)=x2+y2+z2的等值面x2+y2+z2=2，并绘制其梯度场的梯度矢量。

解输入

delta=0.2；

v=-2：delta：2；

［x,y,z］=meshgrid(v)；

u=x.^2.+y.^2+z.^2；［px,py,pz］=gradient(u,delta,delta,delta)；

quiver3(x(1：4：end),y(1：4：end),z(1：4：end),px(1：4：end),py(1：4：end),pz(1：4：end))

p=patch(isosurface(x,y,z,u,2))；

set(p,′FaceColor′,′red′,′EdgeColor′,′none′)；

daspect(［111］)

axistight；

camlight；

lightinggouraud；

输出图形如图7.2所示。图7.2［例3］绘制u(x,y)=ye－x2－y2的等高线及u(x,y)所在梯度场的梯度矢量。

解输入

delta=0.2；

v=-2：delta：2；

［x,y］=meshgrid(v)；

u=y.*exp(-x.^2-y.^2)；

［px,py］=gradient(u,delta,delta)；

［C,h］=contour(x,y,u)；

set(h,′ShowText′,′on′,′TextStep′,get(h,′LevelStep′)*1)

colormapcool

holdon

quiver(x,y,px,py)

holdoff

输出图形如图7.3所示。图7.3

2.散度和旋度计算

［例4］求A(x,y,z)=xy2z2i+z2sinyz3j+x2eyk的散度和

旋度。

解输入

symsxyz

A=［x*y^2*z^2；z^2*sin(y)；x^2*exp(y)］；

DA=jacobian(A,［x,y,z］)结果为

DA=

［y^2*z^2,2*x*y*z^2,2*x*y^2*z］

［0,z^2*cos(y),2*z*sin(y)］

［2*x*exp(y),x^2*exp(y),0］继续输入

pretty(y^2*z^2+z^2*cos(y)+0)

结果为

z2y2+z2cosy

即

divA=z2y2+z2cosy继续输入

pretty(［x^2*exp(y)-2*z*sin(y),2*x*y^2*z-2*x*exp(y),0-2*x*y*z^2］)

结果为

［x2exp(y)－2zsin(y),2xy2z－2xexp(y),－2xyz2］

即

rot

A=(x2ey－2zsiny)i+(2xy2z－2xey)j－2xyz2k［例5］利用MATLAB函数divergence和curl求A(x,y,z)=(x3+z)i+y3j+sin(x+y)k的梯度和旋度，并画出等值面

divA=4，A和rotA的矢量线。解输入

xmin=－7；xmax=7；

ymin=－7；ymax=7；

zmin=－7；zmax=7；

x=linspace(xmin,xmax,30)；

y=linspace(ymin,y