精通matlab科学计算

14编程实现一些偏微分方程的求解。用，而且还能掌握利用对椭圆型偏微分方程、抛物线型偏微分方程和双曲线型14.1偏微分方程概

uA(x,y)2xB(x,y)xyC(x,y)2yf(x,y,u,x,y

剟 xf,y0剟 yu(x,y0)by0(x),u(x,yf)byfu(x0,y)bx0(y),u(xf,y)byf(B24ACB24ACB24AC获得，本章主要讲述它们采用进行的数值求解法。14.2椭圆型偏微分方以下仅以一类特殊的椭圆型偏微分方程Helmhotz方程为例，介绍椭圆型偏微分方程常规Helmholtz2u(x,y)g(x,y)u(x,y)

2u(x,y)

2u(x,

g(x,y)u(x,y)f(x,其自变量取值范围D

剟 xf,y0剟 yu(x,y0)by0(x)，u(x,yf)byfu(x0,y)bx0(y)，u(xf,y)byf(Helmholtzg(xy0Possiong(xy0和f(xy0，则它被称为Laplace为了采用不同的方法求解对区域D进行分割将 xf范围内的x轴等分为M段每段长为：y轴分割成My段每段长为y(yfy0)/My。在Helmholtz方程中，利用三点中心差分估计二阶导数：

2u(x,

|x,y

ui,j12ui,jui,jxjx0jx,yiy02u(x,

|xj,yi

ui1,j2ui,jui1,ui,ju(xj,yi因此，对于自变量区间上的每一个内点(xjyiui,j12ui,jui,j1ui1,j2ui,jui1,jg

i,ji, i,ui,ju(xj,yi),fi,jf(xj,yi),gi,jg(xj,yi在区域D上共有(My1)(Mx1个内点，因此可以建立(My1)(Mx1)个方程。MxMy较大时，联合求解方程组显然不现实，迭代法可以解决这个问题。ui,jry(ui,j1ui,j1)rx(ui1,jui1,j)rxy(gi,jui,jfi,j ui,0bx0(yi),ui,Mbxf(yi),u0,jby0(xi),uM,jbyf(xj

2(x2y2

ry，2(x2y2)rx，2(x2y2)在中编程实现的Helmholtz方程求解程序为：Helmholtz。功能：迭代法求解Helmholtz方程。调用格式：[uxyHelmholtz(fgbx0,bxfby0,byfDMxMyMinErrMaxIter)。其中，f，g为函数名；bx0xx0上的函数值；bxf为在边界xxf上的函数值；by0yy0上的函数值；byf为在边界yyf上的函数值；Mx为沿x轴的区间数；My为沿y轴的区间数；MinErr为误差控制因子；MaxIter[uxy]u(x,y)在(x,y)迭代法求解Helmholtz方程的程序代码如下function[u,x,y]=%解方程：u_xxu_yyg(x,y)u%bx0：在边界x＝x0上的函数值%bxf：在边界x＝xf上的函数值%by0：在边界y＝y0上的函数值%byf：在边界y＝yf上的函数值%D：边界点值，长度为四的列向量，其元素分别为%Mx：沿x轴的区间数%My：沿y轴的区间数%[u,x,y]：方程u(x,y)在(x,y)点的函数值x0=D(1);xf=D(2);y0=D(3);yf=dx=(xf-x0)/Mx;x=x0+[0:Mx]*dx; dy=(yf-y0)/My;y=y0+[0:My]'*dy;Mx1=Mx+1;My1=My+form=u(m,[1Mx1])=[bx0(y(m forn=u([1My1],n)by0(x(nbyf(x(n))];%sum_of_bv=sum(sum([u(2:My,[1Mx1])u([1My1],2:Mx)']));u(2:My,2:Mx)=sum_of_bv/(2*(Mx+My-2));fori=1:Myforj=1:MxF(i,j)=f(x(j),y(i));G(i,j)=dx2=dx*dx;dy2=dy*dy;dxy2=2*(dx2+dy2);rx=dx2/dxy2;ry=dy2/dxy2;rxy=rx*dy2;foritr=1:MaxIterforj=2:Mxfori=2:Myu(i,j)ry*(u(i,j1)+u(i,j1rx*(u(i1,j)+u(i1,j))+rxy*(G(i,j)*u(i,j)-F(i,j)); ifitr1&max(max(abs(uu0MinErr循环结束条件

u0=HelmholtzPossionLaplace14-1】HelmholtzHelmholtz2u(x,y)2u(x,y) 0剟

40剟

4

xu(x,y)x边界u(0yy2u(4y16cos(y)u(x,0x2u(x4x解：可知f(x,y)x2y2，g(x,y) x 中编写函数ex1401来实现求解，代码如下function%调用函数Helmholtz，求满足一定初始条件和边界条件的Helmholtzf=g=x0=0;xf=4;y0=0;yf=4; Mx=30;My=30; bx0=inline('y^2','y'); bxf=inline('4^2*cos(y)','y');by0=byf=D=[x0xfy0yf];MaxIter=400;MinErr=1e-[U,x,y]=Helmholtz(f,g,bx0,bxf,by0,byf,D,Mx,My,MinErr,MaxIter);clf,mesh(x,y,U);14-1U0U43

4 2 0 14-1【例14-1】的Helmholtz方程的数g(xy0Possion方程。此时，只需将程序ex1401.m中的代码“g=inline('sqrt(x)','x','y');”改为“g=inline('0','x','y');”即可。可得其数值解的图形如图14-2所示。若二者均为0，则演变为Laplace方程，即将“g=inline('sqrt(x)','x','y”改为“g=inline('0','x',y');'x','y');，可得其数值解如图14-3所示。U0U43

4 2 0 14-2【例14-1】的Possion方程的数值U0U43

4 2 0 14-3【例14-1】的Laplace方程的数值满足牛顿边值条件的Helmholtzu(x,y)

(Helmholtz方程，将左边界微分方程用差分方程估计，可以ui,1ui,1

(y),

u

(y)x，i1,

My1

ui,0ry(ui,1ui,1)rx(ui1,0ui1,0)rxy(gi,0ui,0fi,0ry(ui,1ui,1

(y)x)r

)

(g

f

i

2ryui,1rx(ui1,0ui1,0)rxy(gi,0ui,0fi,0

yix)MyyMyu0,jry(u0,j1u0,j1)2rxu1,jrxy(g0,ju0,jf0,j

特别地，在边界点(x0y0)u0,02(ryu0,1rxu1,0)rx,y(g0,0u0,0f0,0

(y0)/x

(x0y))在中编程实现满足牛顿边值条件的Helmholtz方程求解程序为：Helmholtz。功能：求解满足牛顿边值条件的Helmholtz方程。调用格式：[u,x,y]=Helmholtz_Newton(f,g,dbx,bx0,bxf,by0,byf,D,Mx,My,其中，f，gdbxbx0xx0上的函数值；bxf为在边界xxf上的函数值；by0yy0上的函数值；byf为在边界yyf上的函数值；D4x0，xf，y0，yf；Mx为沿x轴的区间数；My为沿y轴的区间数；MinErr为误差控制因子；MaxIter[uxy]u(x,y)在(x,y)迭代法求解满足牛顿边值条件的Helmholtz方程的程序代码如下%解方程：u_xxu_yyg(x,y)u%其他条件同Helmholtz函x0=D(1);xf=D(2);y0=D(3);yf=dx=(xf-x0)/Mx;x=x0+[0:Mx]*dx; dy=(yf-y0)/My;y=y0+[0:My]'*dy;Mx1=Mx+1;My1=My+fori=1:Myforj=1:MxF(i,j)=f(x(j),y(i));G(i,j)=fordx2=dx*dx;dy2=dy*dy;dxy2=2*(dx2+dy2);rx=dx2/dxy2;ry=dy2/dxy2;rxy=rx*dy2;form=u(m,[1Mx1])=[bx0(y(m forn=u([1My1],n)by0(x(n sum_of_bv=sum(sum([u(2:My,[1Mx1])u([1My1],2:Mx)']));u(2:My,2:Mx)=sum_of_bv/(2*(Mx+My-2))foritr=form2:(My1- form=2:(Mx1-

forj=2:Mxfori=2:My

u(i,j)ry*(u(i,j1)+u(i,j1rx*(u(i1,j)+u(i1,j))+rxy*(G(i,j)*u(i,j)-F(i,j)); ifitr1&max(max(abs(uu0MinErr循环结束条件u0=14-2】Helmholtz14-1】u|xxy2,u|yx 解： 中编写函数ex1402来实现求解function%调用Helmholtz_Newton解满足一定边值条件的方程：u_xxu_yyf= %构造函数xg=inline('sqrt(x)','x','y'); %构造函数g(x,y)x0=0;xf=4;y0=0;yf=4; Mx=30;My=30; xbx0=inline('y^2','y'); bxf=inline('4^2*cos(y)','y');by0=byf=D=[x0xfy0yf];MaxIter=100;MinErr=1e-clf,mesh(x,y,U)14-4U0U43

4 2 0 14-4牛顿边值的Helmholtz方程数值14-414-1可以发现，两者的区别也主要体现在边界上，这是由于后抛物线偏微分方

2u(x,A

0剟 xf,0剟 u(0,t)b0(t)，u(xf,t)bxf(t)，u(x,0)u0本节的一维抛物线型方程的数值解法有显式前向欧拉法、隐式后向欧拉法将[0xfMxxf/MNtT/N2u(x,A

uk2ukuk

uk1ukA i1 i

uk1

MM

r

在中编程实现的显式前向欧拉法求解抛物线方程的程序为：EF_Euler。调用格式：[uxtEF_Euler(AxfTit0,bx0,bxfMN)。其中，A为抛物线偏微分方程系数；xfxTtit0为在边界t0上的函数值；bx0x0上的函数值；bxfxxfMxNt：方程显示前向欧拉法求解抛物线方程的程序代码如下function[u,x,t]=%解方程Au_xxu_t，0xxf0t%初值u(x,0边界条件u(0,tbx0(tu(xf,t%M：x轴的等分段%N：t轴的等分段dx=xf/M;x=[0:M]'*dx;dt=T/N;t=[0:N]*dt;fori=1:M+1u(i,1)=forn=1:N+u([1M+1],n)=[bx0(t(n));r=A*dt/dx/dx,r1=1-2*r;fork=fori=u(i,k+1)=r*(u(i+1,k)+u(i-1,k))+14-3】显式前向欧拉法求解一维抛物线方程应用实例。求满足下列条件的热传

2u(x,A

u(x,t)，

0)sinxx(,在中编写函数ex1403来实现求解：function%调用函数EF_EulerA= it0=inline('sin(pi*x)','x'); bx0inline('0bxfinline('0');%边界条件xf=2;M=50;T=0.1;N=100;[u1,x,t]=EF_Euler(A,xf,T,it0,bx0,bxf,M,N);[u1,x,t]=EF_Euler(A,xf,T,it0,bx0,bxf,M,N);14-514-6x642U0U2

1 0

x

14-5利用显式前向欧拉法得到的抛物线方程数值解：不稳定区1U0U2

1 0

x

14-6利用显式前向欧拉法得到的抛物线方程数值解：稳定区2u(x,A

uk2ukuk

ukukA i1

(12r)uk

uk

r

若uk，uk的值满足Dirichlet型边界条件，则

(12r)uk

uk1

1 00000000

1u|x0b' ukuk 1b' 等式

(12r)uk

uk1中取k=0 iruk(12r)ukruk i (12r)uk2rukuk12rb

uk12rb'

1

uk 01 101 22 22

1

0

uk

r

uk

M

M

uk

M1 M 在中编程实现的隐式后向欧拉法求解抛物线方程的程序为：IB_Euler。调用格式：[uxtIB_Euler(AxfTit0,bx0,bxfMN)。其中，A为抛物线偏微分方程系数；xfxit0为在边界t0上的函数值；bx0x0上的函数值；bxfxxfMxNt[uxt]u(x,t)在(x,t)隐式后向欧拉法求解抛物线方程的程序代码如下function[u,x,t]=%解方程A1u_xxu_t，0xxf,0t%初值u(x,0边界条件u(0,tbx0(tu(xf,tM：xN：tdx=xf/M;x=[0:M]'*dx;dt=T/N;t=[0:N]*dt;fori=1:M+1,u(i,1)=it0(x(i));forn=1:N+1,u([1M+1],n)=[bx0(t(n));bxf(t(n))];r=A*dt/dx/dx;r2=1+2*r;fori=1:M-1P(i,i) ifi>P(i-1,i)=-r;P(i,i-1)=-fork=2:N+b=[r*u(1,k);zeros(M-3,1);r*u(M+1,k)]+u(2:M,k-1);u(2:M,k)=linsolve(P,b);14-4】隐式后向欧拉法求解一维抛物线方程应用实例。利用隐式后向欧拉法求

2u(x,A

，u0)sinxx(,在中编写函数ex1404来实现求解：function%调用IB_EulerA= it0inline('sin(pi*x)','x');%初始条件bx0inline('0bxfinline('0');%边界条件xf=2;M=80;T=0.1;N=100;[u1,x,t]=IB_Euler(A,xf,T,it0,bx0,bxf,M,N);14-71U0U2

1 0

x

14-7利用隐式后向欧拉法得到的抛物线方程数uk2ukuk

ukukA i1 其左边的差分估计在时间点kt/2的kk-1的中点进行。等式左右的估计点不在同一个时刻。为了消kk+1Auk12uk1uk

uk2uk

uk1uk( i1 i1) ruk12(1r)uk1ruk1ruk2(1r)uk

r

x0/xM2

uk010

2

uk 23 3

2

uk 2

r

M 22rukMuk

k 2

1

u0 2

uk 2

2

uk

3

2

M 22r

(k1)

(k 在中编程实现的Grank-Nicholson法求解抛物线方程的程序为：Grank_调用格式：[uxtGrank_Nicholson(AxfTit0,bx0,bxfMN)。其中，A为抛物线偏微分方程系数；xfxTtit0为在边界t0上的函数值；bx0x0上的函数值；bxfxxfMxNt：方程Grank-Nicholson法求解抛物线方程的程序代码如下function[u,x,t]=%解方程Au_xxu_t，0xxf0t%初值u(x,0边界条件u(0,tbx0(tu(xf,tM：xN：tdx=xf/M;x=dt=T/N;t=fori=1:M+1,u(i,1)=it0(x(i));forn=1:N+1,u([1M+1],n)=[bx0(t(n));bxf(t(n))];endr=A*dt/dx/dxr1=2*(1-r);r2=2*(1+fori=1:M-1P(i,i)=ifi>1P(i-1,i)=-r;P(i,i-1)=-Q(i-1,i)=r;Q(i,i-1)=fork=2:N+u(2:M,k)=linsolve(P,b);14-5】Grank-NicholsonGrank-A0.50剟

Grank-Nicholson14-32u(x,A

u(x,t)，

0)sinxx( 中编写函数ex1405来实现求functionA= it0=inline('sin(pi*x)','x'); bx0inline('0bxfinline('0');%边界条件xf=2;M=25;T=0.1;N=100;[u,x,t]=Grank_Nicholson(A,xf,T,it0,bx0,bxf,M,N);14-8x321U0U2

1 0

0x

14-8利用Grank-Nicholson得到的抛物线方程数值结论：3节中，用三种方法求解一维热传导方程。其中，显式前向欧拉法迭代式的算法思路清晰，容易编程实现，但存在不稳定区域，必须限制r0.5Grank-Nicholson方法利用三对角矩阵进行求解，无不稳定区域，但要求矩阵方程的解，

2u(x,

2u(x,

)u(x,

x xf,y0剟 yf,0剟 u(x0,y,t)bx0(y,t),u(xf,y,t)bxf(u(x,y0,t)by0(x,t),u(x,yf,t)byfu(x,y,0)i0(x,uk

2uk1uk

uk 2uk

uk2ukA(i,j

i, i,j1

i1, i, i1,j)

i, i,ry(uk1uk1)(12r)uk1r

)(12r i, i,j

i,j

i,rx(uk2uk2)(12r)uk2r(uk1`uk1)(12r

i,

i,rxAt

i, i,ryxxfx0,yyfy0,tM

M 在中编程实现的求解二维抛物线方程的程序为：TDE。调用格式：[uxytTDE(ADTixy0,bxytbxfMxMyN)Dx，yTtixy0t0bxyt为边界取值函数；Mxx轴的等分段数；Myy轴的等分段数；N为沿t轴的等分段数；[u,x,y,t]为方程u(x,y,t)在(x,y,t)点的函数值。 function[u,x,y,t]=%初值u(x,y,0边界条件u(x,y,t%Mx，My：x轴和y轴的等分段%N：t轴的等分段dx=(D(2)-D(1))/Mx;x=dy=(D(4)-D(3))/My;y=dt=T/N;t=%初始化forj=1:Mx+1fori=1:My+1u(i,j)=rx=rx1=1+2*rx;rx2=1-2*rx;ry=A*dt/(dy*dy)ry1=1+2*ry;ry2=1-2*ry;forj=1:Mx-1P(j,j)=ry1;ifj>1P(j-1,j)=-ry;P(j,j-1)=-fori=1:My-1Q(i,i)=rx1;ifi>1Q(i-1,i)=-rx;Q(i,i-1)=-fork=u_1=u;t=fori=1:My+1 u(i,1)=feval(bxyt,x(1),y(i),t);u(i,Mx+1)=forj=1:Mx+u(1,j)=u(My+1,j)=feval(bxyt,x(j),y(My+ifmod(k,2)==0fori=2:Myjj=

u(i,jj)=forj=2:Mxii=2:My;+1))+

u(ii,j)=】

A103x0剟

5y0剟

50剟

b(x,y,t)x2sinyy2sinu(xy0)0解：在中编写函数ex1406来实现求解functionex1406()A=1e-3;it0=bxyt=inline('x^2*sin(y)-y^2*sin(x)','x','y','t');D=[0505];T=1000;Mx=50;My=50;N=50;[u,x,y,t]=TDE(A,D,T,it0,bxyt,Mx,My,N);U0U642y

543210 14-9二维抛物线方程的数14.4双曲线型偏微分

2u(x,A

2u(x,t0剟

u(0,t)b0u(x,0)i0

u(xf,t)bxfu(x,0)i' uk2ukuk uk2ukukA i1

uk1

t )2(1r)uk

trAxu11

)(1r)u0i'(x

另，稳定性要求r„1，而算法的精度随着r的增大而提高，因此，选择r1在中编程实现的显式中心差分法求解双曲线型偏微分方程的程序为调用格式：[uxtECD_Wave(AxfTit0,i1t0,bx0,bxfMN)。其中，A为双曲线型偏微分方程系数；xfxTtit0t0i1t0t0bx0x0bxfxxfMxNt[u,x,t]u(x,t)在(x,t)点的函数值。function[u,x,t]=%解方程au_xxu_tt，0<=x<=xf,初始条件u(x,0it0(xu_t(x,0)边界条件u(0,tbx0(tu(xf,t)%M：沿x轴的等分段%N：沿y轴的等分段dx=xf/M;x=[0:M]'*dx;dt=T/N;t=[0:N]*dt;fori=1:M+1u(i,1)=fork=1:N+u([1M+1],k)=[bx0(t(k));r=A*(dt/dx)^2;r1=r/2;r2=2*(1-u(2:M,2)=r1*u(1:M-1,1)+(1-r)*u(2:M,1)+r1*u(3:M+1,1)fork=3:N+-14-7】2u(x,t)2u(x,

0剟

t1，0剟 u0)xx2

u(x,0)u(0,t)0,u(1,t)解： functionex1407()A=1;it0=inline('x-x^2','x');i1t0=inline('0');bx0t=inline('0');bxft=inline('0');xf=1;M=10;T=2;N=[u,x,t]=figure(1), figure(2),clfforn=plot(x,u(:,4*n)),axis([0xf-0.30.3])14-1014-11U0U10y

100

214-10一维波动方程的数0

000-0

0

0

0

000-0

0

0

14-11一维波动方程的二维直观14-114

2u(x,

2u(x,))

u(x,

x xf，y0剟 yf，0剟 u(x0,y,t)bx0(y,t),u(xf,y,t)bxf(u(x,y0,t)by0(x,t),u(x,yf,t)byf

u(x,y,0)i0(x,

u(xy0)i0xy

uk

uk1

ukA(i,j

i, i,j1

i, i1,j)

i, i, i,uk1r

r

r

)uk

i, i,

i,

i, i,rxx2,ry

xxfx0Mx

yyfy0M

tNu11{r

)r )}2(1rr

i'(x,xi, 2 i,j

i,j

i1, i,

在中编程实现的求解二维双曲线型偏微分方程的程序为：Wave2。调用格式：[uxtWave2(ADTit0,i1t0,bxytMxMyN)。其中，A为双曲线型偏微分方程系数；Dx，yTtit0t0i1t0t0bxyt为在边界取值函数；Mxx轴的等分段数；Myy轴的等分段数；N为沿t轴的等分段数；[uxt]u(x,t)在(x,t)点的函数值。function[u,x,y,t]=初始条件u(x,y,0it0(x,yu_t(x,y,0)边界条件u(x,y,t%Mx，My:x轴和y轴的等分段%N：时间轴的等分段dx=(D(2)-D(1))/Mx;x=dy=(D(4)-D(3))/My;y=dt=T/N;t=u=zeros(My+1,Mx+1);ut=zeros(My+1,Mx+1);forj=2:Mxfori=u(i,j)=it0(x(j),y(i));ut(i,j)=adt2=A*dt*dt;rx=adt2/(dx*dx);ry=adt2/(dy*dy);rxy1=1-rx-ry;rxy2=rxy1*2;u_1=u;fork=0:Nt=fori1:My u(i,[1Mx+1])=[bxyt(x(1),y(i),t)bxyt(x(Mx+forj=1:Mx+u([1My+1],j)=[bxyt(x(j),y(1),t);bxyt(x(j),y(My+ifk==fori=2:Myforj=2:Mxu(i,j)=0.5*(rx*(u_1(i,j-1)+u_1(i,j+fori=2:Myforj=2:Mxu(i,j)=rx*(u_1(i,j-1)+u_1(i,j+1))+ry*(u_1(i-1,j)u_1(i+1,j))+rxy2*u(i,j)-

u_2=u_1;u_1=14-8】

2u(x,

2u(x,))

u(x,0剟

20剟

u(0,y,t)0,u(2,y,t)u(x,0,t)0,u(x,yf,t)u(x,y,0)cos(x)sin(y),u(x,y,0)解： functionit0=i1t0=inline('0','x','y');bxyt=inline('0','x','y','t');A=0.5;D=[0202];T=2;Mx=40;My=40;N=40;[u,x,y,t]=Wave2(A,D,T,it0,i1t0,bxyt,Mx,My,N);14-121U0U2

1 0

1x

214-12二维波动方程的数14.5有限元有限元法是求解偏微分方程数值解的另法，它既能处理形如[X,Y,Z]的规则边2u(x,y)2u(x,y)g(x,y)u(x,y)

f(x,Du(x,y)b(x,DNs个区域，{S1

S}NsN在边界上确定Nn

{N1,

Nb}，{Nb1,

Nn}n(x,y){n,s,s1, ,Ns}(x,y)n,s(x,y)pn,s(1)pn,s(2)xpn,s(3)这样的定义可以保证nn10。偏微分方程的解可以由n线性表示： 1 2

u(x,y)cT(x,y)

cnn cnn cnncT nNb1

N]T,c1[c1,

cN2[N1,N

N]T,c2[cN1,cN

cN s(x,y)cnn,s(x,y)cn(pns(1)pns(2)xpns(3) 根据边界条件确定c1 2,s1,s 2,s1,s

2,s

1,s

2,s

1,s

g(x,y)T 2,s2,s 2,s2,s

2,s

2,s

2,s

2,s

g(x,y)T

dA1c1f(xs,ys

A1c1A2c2 f(xs,ys)2,sSsI

u(x,y)

u(x,y)

g(x,y)u2(x,y)2g(x,y)u2(x,y)

cTTccTTcg(x,y)2Tcf(x,y)2R 首先通过函数fem_basis_ftn产生基函数，使得nn10。在中编程实现的产生基函数的程序为：fem_basis_ftn。调用格式：pfem_basis_ftn(N,S)。其中，N为节点坐标矩阵；S为子区域的节点列表矩阵，第i行为第i个子区域内的所有节点P为基函数ftnphi_i系数。functionp=%p(i,s,1:3基函数ftnphi_i系数，共s%N(n,1:2)第n个节点的x和y坐%S(s,1:3)第s个子区域的节点#sN_n=size(N,1);%总共节点数N_s=size(S,1);%总子区域个数forn=1:N_nfors=1:N_sfori=A(i,1:3)=[1b(i)S(s,i)n第n个基函数仅在节点n

fori=1:3,p(n,s,i)=pnt(i);通过函数fem_coef为每一个节点和子区域构造基函数n,s(xy)，并求得系数向量c；在中编程实现产生的为每一个节点和子区域构造的基函数为：fem_coef。调用格式：[U,cfem_coef(fgpcNSN_i)。其中，f为拉斯方程等式右边的函数f(x,y)；g为拉斯方程等式左边方程g(x,y)；c为边界点的取值；NS为子区域的节点列表矩阵，第i行为第i个子区域内的所有节点N_ifunction[U,c]=%p(i,s,1:3基函数ftnphi_i系%c.1100边界节点取1，内节点取%N(n,1:2)第n个节点的x和y坐%S(s,1:3)第s个子区域的节点%N_i%U(s,1:3)每个区域的p1p2(s)xp3(s)y的坐标N_n=size(N,1);%总共节点数N_ssize(S,1总共子区域N_b=N_n-N_i;fori=forn=1:N_nfors=1:N_sxyN(S(s,1N(S(s,2N(S(s,3),:))/3重p_vctr=[p([in],s,1)p([in],s,2)p([i[1xy]'*p_vctr(2,:)*[1xy]';dS(s)=det([N(S(s,1),:)1;N(S(s,2),:)1;N(S(s,3),:)ifn==1,tmpf(s)=-f(xy(1),xy(2))*p_vctr(1,:)*[1xy]';A12(i-N_b,n)=d(i-N_b)=d=d-A12(1:N_i,1:N_b)*c(1:N_b)';c(N_b+1:N_n)=A12(1:N_i,N_b+1:N_n)\d;fors=1:N_sforj=1:3,U(s,j)=c*p(:,s,j); function%N11;11;11;-11;0.2 N10;01;12;21;1.2 N_n=size(N,1); S=[125;235;345;14 N_s= figure(1),clffors=1:N_snodes=[S(s,:)S(s,1)];fori=1:3holdonp=%x01;xf1;y01yf x00;xf2;y00;yf figure(2),clfMx=50;My=50;c=[0123 %dx=(xf-x0)/Mx;dy=(yf-xi=x0+[0:Mx]*dx;yi=y0+i_ns[1234 foritr=i_n=i_ns(itr);ifitr==1fori=1:length(xi)forj=1:length(yi)Z(j,i)=fors=ifinpolygon(xi(i),yi(j),N(S(s,:),1),N(S(s,:),2))>

subplot(321mesh(xi,yi,Z);title(itr)节点1c1=zeros(size(c));c1(i_n)=1;subplot(320+itr) %节点2~5的基函c=[012314-1314-1411110000 0 14-13节点和子 1110214-14基函数及其线性组图【例14-9】有限元法求解偏微分方程应用实例。采用有限元法求解如下拉斯方2u(x,y)2u(x,y)

1剟 1，1剟

f(x,11,(x,y)(0.5,ff(x,y)1,(x,

(0.5,u(xy)0

0,(x,y)注意到只有在内点(0.50.5),(0.50.5附近，函数f(xy的值不为零，其余地方均为零。u(xy的值会存在剧烈的变化：出现波峰或者波谷。functionN=[-10;-1-1;-1/2-1;0-1;1/2-1;1-1;10;11;1/21;0-1/21;-11;-1/2-1/4;-5/8-7/16;-3/4-5/8;-1/2--1/4-5/8;-3/8-7/16;00;1/21/4;5/87/16;3/45/8;1/25/8;1/45/8;3/87/16;-9/16-17/32;-7/16-17/32;-1/2-7/16;9/1617/32;7/1617/32;1/2N_b= S=[11112;11119;101119;4519;5719;567;1215;2331517;3417;41719;131719;11319;11315;7822;8992224;91024;101924;192024;71920;72022;1314141516;161718;202125;212223;232425;1426162627;182728;212931;232930;25302627282930 f11='-(norm([xy]+[0.50.5])<0.01)-(norm([xy]-[0.50.5])<N_n= N_i=N_n- c= %边界节点和节点的取p=[U,c]=figure(1),clf,trimesh(S,N(:,1),N(:,2),c)N_s=size(S,1); x0=-1;xf=1;y0=-1;yf=1;Mx=16;dx=(xf-x0)/Mx;xi=x0+[0:Mx]*dx;My=16;dy=(yf-y0)/My;yi=y0+[0:My]*dy;fori=1:length(xi)forj=for