高精度公选课(奇 数)

一．二．答：工程中常用的数值分析方法：1.有限元法（FEM）优点：对于不规则几何区域的适应性好。

为求解大型物体的计算带来了方便，并对问题给出了宏观的把握，从而更好地提出解决措施。

缺点：需要采用高阶单元或采用较小尺寸的单元，极大降低分析问题的效率。

2.有限差分法（FDM）优点：是求得偏微分方程数值解的最古老方法，对于简单的几何形状中流动与传热问题也是一种最容易实施的方法。

在规则区域的结构化网格上，有限差分法是十分简单而有效的，而且很容易引入对流项的高阶格式。

缺点：离散方程的守恒特性难以保证，而最严重的缺点则是对不规则区域的适应性差。

3.边界元法（BEM）优点：可以使求解问题的空间维数降低一阶，从而使计算工作量及其所需计算机容量大大减小。

缺点：需要已知所求解偏微分方程的格林函数基本解，虽然对不少偏微分方程这个基本解已找出，但对Navier-Stokes方程这个非线性偏微分方程，至今尚未找到基本解。

4.有限分析法（FAM）优点：有限分析法是80年代初发展起来的一种数值方法，它可以克服高Reynolds 数下有限差分法及有限容积法的数值解容易发散或振荡的缺点。

缺点：其计算工作量较大，对计算区域几何形状的适应性也较差。

5.有限容积法(FVM）优点：用有限容积法导出的离散方程可以保证具有守恒性（只要界面上的插值方法对位于界面两侧的控制容积式一样的即可），对区域形状的适应性也比有限差分法要好。

是目前应用最普遍的一种数值方法。

缺点：有限体积法在寻求控制体积的积分时，必须假定值在网格点之间的分布。

6.谱分析法（Spectral method，SM）优点：可以获得高精度的解，目前只是在比较简单的流动与传热问题中已经应用的比较成功。

缺点：不适宜用于编制通用程序。

7.数值积分变换法（Integral transformation method，ITM）优点：计算精度可以较高，降维问题是一个常微分方程，有成熟的数值方法可以采用。

缺点：这种方法不容易形成通用程序，特征值问题的选取有一定的任意性，对非线性强烈的问题，计算工作量比较大。

三．答：1.重心插值配点法的基本计算步骤：首先，利用重心插值微分矩阵可以得到微分方程和边界条件的离散形式；其次，施加边界条件：边界条件施加方法主要有消去法、置换法和附加法；再次，对于高阶边值问题，需要对边界条件做一些特殊的处理，以保证计算结果的精度和可靠度。

2.边界条件的离散：对于两点边值问题，微分方程为：()+=+),(''(3.1)<)'(qxaxbfypxyy<x边界条件可表示为B b y B b y B A a y A a y A =+=+)(')(,)(')(2121 (3.2)式中：)()()(x f x q x p 、、均为区间[a,b]上的连续函数；B B B A A A 、、、、、2121均为常数。

将区间[a,b]离散为<<=21x x a ···b x n =<，令,,21y y ···n y ，为未知函数)(x y 在离散节点,,21x x ···n x ，上的函数值，利用重心型插值近似未知函数)(x y 为∑==nj j j y x L x y 1)()( (3.3)利用插值公式(3.3)和微分矩阵的记号,边界条件离散为B y D B y B A y D A y A i nj nj n j n j j=+=+∑∑==1)1(211)1(1211, (3.4) 3.边界条件的施加方法：微分方程(3.1)的离散形式为f y QI PD D =++][)1()2(简记为：f Ly = (3.5)边界条件施加方法主要有以下三种：(1).消去法消去法施加边界条件的方法是，从边界条件离散方程中解出n y y 、1，也就是将n y y 、1用,2y ···1-n y ，表示。

将n y y 、1代入离散方程组(3.5)的第二至第n-1个方程，得到关于,2y ···1-n y ，的方程组求解得到两点边值问题的数值解。

(2).置换法置换法施加边界条件的方法是，利用边界条件离散方程(3.4)的两个方程分别置换方程组(3.5)第1，第n 个方程，实现边界条件的施加。

(3).附加法附加法施加边界条件的方法是，将边界条件离散方程(3.4)的两个方程附加到方程组(3.5)中，构成n 个未知量n+2个方程的方程组求解，实现边界条件的施加。

四．P1.clcclearn=15;a=0;b=1;x=-cos((0:n)*pi/n);x=x*(b-a)/2+(b+a)/2;x=x';D=barylag(x,3); D2=D(:,:,2);D1=D(:,:,1);D3=D(:,:,3);I=eye(size(D1));n=n+1;L=D3-diag(x);f=(x.^3-2*x.^2-5*x-3).*exp(x);ye=x.*(1-x).*exp(x);L(1,:)=I(1,:);f(1)=0;L(n,:)=I(n,:);f(n)=0;L(2,:)=D1(1,:);f(2)=1;yy1=L\f;err1=norm(yy1-ye);rerr1=norm(yy1-ye)/norm(ye);L=D3-diag(x);f=(x.^3-2*x.^2-5*x-3).*exp(x);L(n+1,:)=I(1,:);f(n+1)=0;L(n+2,:)=I(n,:);f(n+2)=0;L(n+3,:)=D1(1,:);f(n+3)=1;yy2=L\f;err2=norm(yy2-ye);rerr2=norm(yy2-ye)/norm(ye);err=[err1 err2 ]rerr=[rerr1 rerr2 ]plot(x,ye,'-ok',x,yy2,'-^b')err =1.0e-13 *0.6923 0.1765rerr =1.0e-13 *0.6796 0.1732 00.10.20.30.40.50.60.70.80.91-0.0500.050.10.150.20.250.30.350.40.45图1.数值解与解析解图形00.10.20.30.40.50.60.70.80.91-5-4-3-2-112x 10-14图2.绝对误差分布图P2.clcclearN=15;T=4;x=-cos((0:N)*pi/N);x=T*(x+1)/2;D=barylag(x,4); D2=D(:,:,2);D1=D(:,:,1);I=eye(size(D1));L=D2+D1-2*I;f=(2*x)';L(N+1,:)=I(1,:);f(N+1)=0;L(N+2,:)=D1(1,:);f(N+2)=1;y=L\fyy=exp(x)-0.5*exp(-2*x)-x-0.5;err=max(abs(y-yy'))rerr=max(abs(y(2:N+1)-yy(2:N+1)')./y(2:N+1))plot(x,yy,'-k',x,y,'.k')y1=D1*y;ye1=exp(x')+exp(-2*x')-1;errd1=max(abs(y1-ye1))rerr1=max(abs(y1(2:N+1)-ye1(2:N+1))./ye1(2:N+1))figureplot(x,ye1,'-k',x,y1,'.k')y2=D2*y;ye2=exp(x')-2*exp(-2*x');errd2=max(abs(y1-ye1))rerr2=max(abs(y2(2:N+1)-ye2(2:N+1))./ye2(2:N+1))figureplot(x,ye2,'-k',x,y2,'.k')err =4.3875e-09rerr =5.8501e-10errd1 =5.7623e-09rerr1 =7.2107e-10errd2 =5.7623e-09rerr2 =2.6458e-0800.51 1.52 2.53 3.54-100102030405060图3.数值解与解析解图形 00.51 1.52 2.53 3.54-2024681012141618x 10-8图4.绝对误差分布图P3.clc;clearm=15;n=m;br=5;a=-1;b=1;c=0;d=1;x=-cos(pi*(0:m)'/m);x=(b+a)/2+(b-a)*x/2;y=-cos(pi*(0:n)'/n);y=(c+d)/2+(d-c)*y/2;%x=linspace(-0.5,0.5,m+1);x=x';y=linspace(-0.5,0.5,n+1);y=y'; m=length(x);n=length(y);Im=eye(m);In=eye(n);I=eye(m*n); C=barylag(x,2);C1=C(:,:,1);C2=C(:,:,2);D=barylag(y,2);D1=D(:,:,1);D2=D(:,:,2);[X,Y]=meshgrid(x,y);X=X(:);Y=Y(:);L=kron(C2,In)+kron(Im,D2);F=8*X;u=inline('x.*(x.^2+y.^2)+exp(3*x).*cos(3*y)','x','y');bcxl=[1:n];bcxr=[(m-1)*n+1:m*n];bcyd=[1:n:m*n];bcyu=[n:n:m*n];L(bcxl,:)=I(bcxl,:);F(bcxl)=u(a,y);L(bcxr,:)=I(bcxr,:);F(bcxr)=u(b,y);L(bcyd,:)=I(bcyd,:);F(bcyd)=u(x,c);L(bcyu,:)=I(bcyu,:);F(bcyu)=u(x,d);Ul=L\F;Ue=u(X,Y);Ue=Ue(:);Er=norm(Ul-Ue);Rer=Er/norm(Ue);[Er Rer]ul=reshape(Ul,n,m);ue=reshape(Ue,n,m);[x,y]=meshgrid(x,y); surf(x,y,ue); figure; surf(x,y,ul-ue)ans =1.0e-09 *0.1773 0.0016 -1-0.50.5100.51-20-10102030图5.解析解函数图形 -1-0.50.5100.51-4-2246x 10-11图6.误差分布图P4.clear;clcc=1;a=0;b=1;l=1;Tf=4;m=15;n=m;dx=8;dt=10;x=-cos((0:m)*pi/m);x=x';x=0.5*(b-a)+0.5*(b-a)*x;t=-cos((0:n)*pi/n);t=t';t=0.5*(t+1)*Tf;%x=linspace(a,b,m+1);x=x';t=linspace(0,Tf,n+1);t=t';m=length(x);n=length(t);Im=eye(m);In=eye(n);I=eye(m*n); Dx=barylag(x,2);Dx1=Dx(:,:,1);Dx2=Dx(:,:,2);Dt=barylag(t,2); Dt1=Dt(:,:,1);Dt2=Dt(:,:,2);Dx=baryrat(x,dx,2);Dx1=Dx(:,:,1);Dx2=Dx(:,:,2);Dt=baryrat(t,dt,2); Dt1=Dt(:,:,1);Dt2=Dt(:,:,2);[X,T]=meshgrid(x,t);X=X(:);T=T(:);L=-0.5*diag(X.^2)*kron(Dx2,In)+kron(Im,Dt1);F=zeros(m*n,1); bcxl=[1:n];bcxr=[m*n-n+1:m*n];bct0=[1:n:m*n];L(bcxl,:)=I(bcxl,:);F(bcxl)=0;L(bcxr,:)=I(bcxr,:);F(bcxr)=exp(t);L(bct0,:)=I(bct0,:);F(bct0)=x.^2;Uc=L\F;Ue=X.^2.*exp(T);Er=norm(Uc-Ue);Rer=Er/norm(Ue);[Er Rer]uc=reshape(Uc,n,m);ue=reshape(Ue,n,m);[x,t]=meshgrid(x,t); surf(x,t,ue); figure; surf(x,t,uc-ue)ans =1.0e-07 *0.1239 0.0006 00.20.40.60.8101234102030405060图7.解析解函数图形00.20.40.60.8101234-3-2-11234x 10-9图8.绝对误差分布图。