某大学计算方法课件-第七章常微分方程边值问题数值解法.pdf

《某大学计算方法课件-第七章常微分方程边值问题数值解法.pdf》由会员分享，可在线阅读，更多相关《某大学计算方法课件-第七章常微分方程边值问题数值解法.pdf（63页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第七章常微分方程边值问题数值解法一个由常微分方程及其边界条件构成的定解问题称为常微分方程边值问题，其常常出现在流体力学、材料力学、波动力学、电磁学及现代控制理论等学科中。一般而言，常微分方程边值问题的解析求解是非常困难的。为了解决这类实际问题，人们通常采用数值解法。本节将主要介绍71维常微分方程边值问题f，（力）=/（力。（力），i C（Q 2）,的打靶法、有限差分法及有限元法。7.1打靶法从形式上看边值问题与初值问题之间的差别仅在于定解条件的不同，毫无疑问这两类问题有着紧密的关联。事实上，基于同一微分方程的边值问题与初值问题有着共同的通解形式。因此，借助于初值问题数值算法来构造边值问题的数值

2、算法是一个很直接的思路，而前者有很多有效的数值算法可供选择。打靶法的基本思想是从不同的初值出发用已知的初值问题数值算法解相同的微分方程直到解的轨迹“打”到给定的右边界点。为简单起见，本节我们仅考虑如下几维线性边值问题f yr=A(t)y+q(t),t G (Q,6),BaM +Bby(b)=7 7,()其中4 以和场为X T Z的矩阵q和为几维向量，而所用方法也可拓展到非线性边值问题。7.1.1简单打靶法对于一个常微分方程而言，不同的初值条件对应于不同的解的轨线，即解是依赖于初始条件的。令沙依。)为 方 程/=力,)满足初值条件g(Q)=。的解，若g(力；。)同时也满足边值问题(2)的边界条件

3、，那么g(1;c)正好就是边值问题(2)的解。简单打靶法要做的工作正是如何寻求这样的初值c,使得BQC+Bby(b c)=q(3)根据常微分方程理论，我们知道线性边值问题的解可表示为y(t)=Y(t)c+v(t),t e(a,6),(4)其中Y )，。分别为基解矩阵及特解，为一参数向量，且满足Y，=4 K 力 G (Q,b),Y(Q)=/;(5)v =A(t)v+q，O(Q)=a,(6)这里/为级单位矩阵，a G史通常可取为零向量。通过计算初值问题(5)-(6),我们可以得到边值问题的解的表达式(4)。为了确定参数向量c,将(4)代入问题的边界条件可得Q c=力，(7)其中Q ：=Ba+B6y

4、(5),r j ：=r j-Bav(a)-Bhv(b).若边值问题(2)存在唯一解，则Q 是非奇异的，因此c 是存在唯一的。当a 为零向量时，。正好为边值问题的初值条件，在程序中算出c 后，我们再次计算一次初值问题得到原边值问题的数值解。_简单打靶法的计算程序如下，其中bas说上八和par税 分 别 表示函数4 1汝和4“g+q(,B a,6b与况a表示边鼻参数，t s p a n 为求解区间。算法7.1 简单打靶法_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _t,Y=shooting(basisfun,partifun,Ba

5、,Bb,tspan,e ta)n=Iength(Ba);I=eye(n);T=zeros(n);t,y =ode45(partifun,tspan,zeros(n,1 )；e ta 1 =e ta Bb*y(length(t),:)for j=1:nt,y=ode45(basisfun,tspan,1(:,j);T(:,j)=y(length(t),:)5;endQ=Ba+Bb*T;c=Qe ta 1 ;t,Y=ode 45(partifun,tspan,c);例7.1试分别取X=1,7,15,3 0,应用简单打靶法求解 0,1上的边值问题/0 1 0/00 0 1 什 02A3 A2 2A/

6、(7r3+A27r)s i n(7rZ)+(2A3+2A7r2)COS(TT)/I 0 0/0 0 0/(e-2A+2e-A+3)/(e-A+2)0 0 0 7/(0)+1 0 0 7/(1)=00 0 0/0 1 0/(A(3 -e-A)/(e-A+2)要 求 画 出 数 值 解 和 精 确 解 的 对 比 图，方程精确解为g=（，”,（力）丁，其中=g A(i-l)+e2A(i 1)+ext2+e-x+COS(7T力).解 分 别 取1,7,1 5,3 0,直接应用简单打靶法7.1可计算出相应的数值解。Z521.510.50-0.500.5 1t图 7.1a A =1.21.51n 0.5

7、0-0.5-10 0.2 0.4 0.6 0.8 1t图 7.1b A=7.20151050-50 0.2 0.4 0.6 0.8 1t图 7.1c A =15.43210.八20 x 105n0 0.2 0.4 0.6 0.8 1t图 7.1d A =30.从上图可见，随着M直的增大，打靶法的误差也增大，出现该现象的原因在于：一方面，相应的矩阵Q=&+石 万 越 来 越 病态；另一方面，直愈大，相应初值问题的刚性也愈强，从而导致数值解在人取相对大的值时会出现较大误差。7.1.2多重打靶法如果初值问题的积分区间较长，简单打靶法的计算效果将会降低。为克服简单打靶法的这一缺点，下面我们介绍多重打靶

8、法，其基本思想是首先将求解区间作网格分划：a=力 N_i =A然后在每个小区间比T 句上构造微分方程的逼近解，最后由解的连续性及边界条件可将这些数值解修补在一起得到一个全局解。考虑线性问题(2),在区间期_ 1,切上，边值问题的解可以表示为y(t;Ci-i)=匕Q 1 +4，。=1,2,，N,(8)其中匕(力)，3分别为基解矩阵及特解，其满足Y：=4 1)匕，匕(力I)=I,(9)3 9 i)=0.(10)计 算 初 值 问 题(9)-(10)得区间比_i,句上边值问题的解的表达式(8)。为了确定参数向量c：=(%c f,，C(1)T，将(8)代入问题的边界条件得g 匕(幻4一1=%(切，1

9、Z?/-1,(11)Bac o +BIYNCN-X=T)B i)V N(b).(12)(11)-(12)可写成如下线性系统Lc =r,(13)其中/-H(力 i)I B(力 2)/L=.,.,YN-ID N-I)I BQ BIYN Jr =(。1(右)2(力2),，n -B wN(b).解线性系统(13)即得网格点力上的数值解纳=Q。多重打靶法求解线性边值问题的计算程序如下，其中输入量表示意义与简单打靶法基本相同，差别仅在力帆es/i为求解区间上的全体网格点。算法7.2 多重打靶法_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

10、Y=multishooting(basisfun,parti fu n,Ba,Bb,tmesh,eta)n=length(Ba);I=eye(n);T=zeros(n);N=length(tmesh)1;L=zeros(N*n);r=zeros(N*n,1);Y=zeros(N+1,n);for k=1 :Nt,y J=ode45(partifun,tmesh(k),tmesh(k+l),zeros(n,1);r(n*(k l)+l:n*i)=y(length(t)for j=1:nt,y J=ode45(basisfun,tmesh(k),tmesh(k+1),1(:,j);T(:,j)=y

11、(length(tendL(n*(kl)+l:n*k,n*(k 1)+1:n*k)=T;if k小L(n*(k l)+l:n*k,n*k+l:n*(k+l)=I;endendL(n*(N-l)+l:n*N,1 :n)=Ba;L(n*(N-1)+1:n*N,n*(N1)+1:n*N)=Bb*T;v=r(n*(N1)+1:n*N);r(n*(N-1 )+1 :n*N)=eta Bb*v;c=Lr;for i=1 :NY(i,:)=c(n*(i l)+l:n*i)endY(N+1,:)=(T*c(n*(Nl)+l:n*N)+v)将多重打靶法应用于例7.1中1=30的情形，可计算出相应的数值 解（见图

12、7.2）。由图7.2可见，在计算精度方面，多重打靶法比简单打靶法具有明显的优势。21n 0.50-0.51.5-10 0.2 0.4 0.6 0.8 1t图7.2例7.1。=30）的精确解（虚线）及其多重打靶解（点戈IJ线）.7.2有 限 差 分 法有限差分法也是一种基于初值问题数值方法的边值问题方法，但其不同于打靶法，有限差分法在利用初值问题数值方法离散边值问题后，将求解区间上所有网格点处的数值解作为离散格式的未知量，从而通过求解表征该格式的线性或非线性代数方程组而获得边值问题的数值解。在本节，我们仍然考虑一阶微分方程系统的边值问题(1)-(2)。为了获得边值问题的系列离散格式，我们对求 解

13、 区 间 作 如 下 剖 分：n ：Q=力1 1 2力N _ 1 且记步长h i =t j 纳为g(切的逼近值。7.2.1中点法与梯形法在本节，我们先考虑简单的有限差分格式，即中点法和梯形法。分别应用中点法和梯形法到非线性边值问题（1）-（2）,我们依次可得其问题的离散格式yt yi-i(友 一1/2,1(%+Vi-1，Q =1,2,.,N),g(y0,yN)=0(14)产t =如+/佐1,统 一1）（分=1,2,，N）,。（如,yN）=0.（15）上 述 计 算 格 式 可 视 为 以 向 量J=（端，端，弧）T为未知量的（N +1）维非线性方程组，其可用Ne wt o n迭代法求解。为叙述

14、简单见，这里我们仅给出梯形法的Ne wt o n迭代格式，中点法的Ne wt o n迭代格式类似可得。记Nh Vi：=纳 统）+/侑1,%1）,i=l,2,，NF：=（Mg即式,N h曲g g 加丁）丁,则中点法的Ne wt o n迭代格式为卜 旷）E=-F 铲）,r+1=r +,m =0,l,2,.-,M（16）其中尸（。为歹。的Ja c o b i 矩阵，加为迭代次数。若引入记号：e=（品 ，册,4 切=力2 n&=&（鳄,第）,Bb=次 端 第）则迭代格式（16）可进一步写成（幻金 +4 1 一1后一1 =N hg.N tol)for j=1 :Nhj=tmesh(j+1)tmesh(j

15、);A=feval(prfun,tmesh(j)+hj/2,(y(j,:)+y(j+l,:)/2);L(n*(j l)+l:n*j,n*(j 1)+1:n*j)=1/hj*1 A/2;L(n*(j l)+l:n*j,n*j+l:n*(j+1 )=1 /hj*1 A/2;r(n*(j-l)+l:n*j)=(y(j+1 y(j,:)/hj+feval(rfun,tmesh(j)+hj/2,(y(j,:)+y(j+1 ,:)/2);endL(n*N+l:n*(N+l),l:n)=feval(pgfu nl,y(l,:)5,y(N+l,:)5);L(n*N+l:n*(N+l),n*N+l:n*(N+l

16、)=feval(pgfun2,y(1,:)5,y(N+lr(n*N+l:n*(N+1)=feval(gfun,y(1,:)5,y(N+1 ,:)*);xi=Lr;for k=1 :N+1y(k,:)=y(k,:)+x i(n*(k l)+l:n*k);endend end上 述 程 序 中，输 入量rfu n和p r f u n分 别 表 示 函 数g f un，p g f un l，p g f un 2 分别表示函数 q(%o),呼 2迦 弗t m e s h 表示求解区间上所有网格点同工1,.工N，y O表示网格点上的初始迭代值。例7.2试取步长九=1/(10*2Z)(i=0,1,2,3)

17、,分别应用中点法和梯形法求解非线性两点边值问题力oQ,_一O,),康e?-=+)/-/O(1X要求给出数值解的全局误差及收敛阶，其方程的精确解为解 令u(t)2 Incosh(11/2)6*/2cosh(8/4),0 =1.5171646.则其边值问题可转化为一阶形式yr=O 力 IB a y 6+牖 =0.取步长九=1/（10*2，）（，=0,1,2,3）,初始迭代值为零，分别应用中点法和梯形法求解该边值问题可计算出其数值解的全局误差及收敛阶（见表7.1）,其中全局误差及收敛阶的表达式依次为/八 、I 1 err(h)err(ri):=max (七 一%,p logo-J oi；v1 k 小

18、 82 167 Td/2)表7中点法和梯形法解例7.2的数值结果.h中点法梯形法err(h)Perr(h)P1/10 6.3341e-4-5.8442e-4-1/201.5938e-41.99071.4632e-41.99791/40 3.9867e-51.99923.6594e-51.99941/80 9.9698e-61.9996 9.1516e-61.9995由表7.1可见，中点法和梯形法均只有2阶精度，并且二者的计算效果相近。7.2.2 Runge-Kutta方法由上可知，中点法和梯形法是二种低精度方法。为获得高精度的边值问题方法，我们考虑修改求解初值问题的Runge-Kutta方法，

19、使其适再于边值问题-（2）。将Runge-Kutta方法应用于边值问题可得Nhfi3：=/o-f=0,l i N,1 j s,SNhVi：=%Vi-i-hi bjfij=0,1 i N g(y&yN)=0,)=1(18)这里玲=+吗。该方程组中的未知向量y：=（诵J f,对，届 ，不，我）fi：=（总 虑 ，忌）,可通过Ne wton迭代法获得，其计算格式为F (=F句 ym+1=r+,m=(19)其中尸为9 3的Jacobi矩阵，机为迭代次数，且尸：=（N/J。Nh V，N h f Nh 状、,Nh代,Nh V g 0,VNY、人 r r-i r r r r r r iNdi=（Nhf l

20、N麻,，明博）丁.若进一步引入记号：=（品 源 卓 ，苏8）7=（4，脸T,Q =-此人sB a =g y：y那,为=心 以,或），Ag）=弘 熄 +儿 口,1=1（4（力 ）（。114/讥），%=：，用=f ：T4（友 s）J QS1A（力 谢）则迭代格式（19）可写成 Q1S4(力 讥)、:,。=(瓦/,C LssA.(tis I(%=Q?:,1 i N,&-1 h,Da,=0,1 S，S N)&a+&切=-g(鳄,。口 一+1 =y+&-0 2 T V,，承+1=印 +%I i T V.(20)既然(20)中的第一个方程可写为口=叱一十%1+叱一崂力1 2 N,因此将其代入(20)中的第

21、二个方程可得&-1=陌，i i N,其中八 口=1+h i D W V i,n =h i D W Q i.从而迭代格式(20)可写成&一-1=片，1 z A,rBa&+1+B b&N=-g(y/,y用,=r+6：,o z 7 v,)举+1=优+叱-1%1 +用一七力l i tol)for i=1 :Nh _i=tme sh(i+1)tme sh(i);for j=1:st _ i j=tme sh(j)+h _ i*rh o(j);y_ij=y(i+h _i*kron(alpha(j,:),I)*f(s*n*(il)+l:s*n*i);A _ij=fe val(prfun,t_ij,y _i

22、j);V_i(n*(j-1)+1:n*j,:)=A_ij;W_i(n*(j-1)+1:n*j,:)=kron(alpha(j,:),A_ij);Q_ij=feval(rfun,t_ij,y _ij)f(s*n*(i l)+n*(j 1)+.1:s*n*(i l)+n*j);Q_i(n*(j-l)+l:n*j,l)=Q _ij;endwu e y e(s*n)h1i*w i)uuw(1)*V；p TW i(l)*Q i;U(二i)H U i;p(：i;r-(n*(i l)+l:n*i。n*(il)+=n*ill(I+hL.*k r o n(b e f a。I)*w i(l)*v i)L(n*(i

23、 L)+l:n*i.n.i+lIn A i+cvI-r oA i 1)+1-n*ihli*k r o n(b e f a,I)*w i(l)*Q i;e n dL(n*N+l:n*(N+l)，一:n)u f e v a l(p g f u nl“y(一：)，y(N+l“：)”)；L(n*N+l-n*z+l)，n*N+l:n*z+l)=f e v a l(p g f u n5,y(l，y(N+l:.)；r(n*N+1n*z+l)nf e v a l(g f u n“y(l,:)r yz+lJ)xiuL-r”f o r k u l:NyQ lTy.:)+x i(n*(k l l)+=n*k)二f(

24、s*n*(k l l)+l-s*n*k n f(s*n*(k l l)+l:s*n*k)+:.u(r:,k)*x i(n*(k l l)+=n*k)+p(：.k)；e ndy(N+l.)u y(N+l.)+x i(n*N+l j n*(N+l);e nde nd分别取步长九=1/(10*2，)(，=0,1,2,3),应用二级RadauIA方法和RadauIIA方法于例7.2,可获得其数值解。在表7.2中，我们给出了数值解的整体误差e rr=maxQi:f (i+j,Vi+j),i=ki,.,N k2,(22)j=-ki j=ki其中 A；1+k2=k,to=a,ti=t i-1 +h,h=(b

25、-a)/N,y,为 y g的逼近。为 了 求 解 方 程 组(2 2),我 们 仍 需 要k个 逼 近值班加，1，以-%+1，以，其可由k-1个与(22)同阶的线性多步公式ki kif y i+j =hf吃J gj,%+/),，=1,岛 一 1,(23)j=i j=ik ikii+jl/N-k+i+j=h f (tN-k+i+j l/N-k+i+j，=N 一 卷+1,.,Nj=T j=-i(24)及边界条件g o,U N)=O(25)给出。至此，我们即得到求解边值问题(1)-(2)的线性多步法(22)-(25)。为实现格式(22)-(25)的计算，我们将(22)-(24)写成矩阵形式其中y=(

26、湍,对，,成 咒 F=(f也,y o F fg 班产,V N)T)T,Qg QIQ。a。(s-fc2+l)Q l -Qka尸2+1).a f f+i)矩 阵B与4 有 相 同 形 式，仅 需 将 A 中 的 内,a?分别代之以如 邛即可获得。联 立(26)与(25),并用Ne wt o n迭代法解之，即可获得彷侑问题(1)-(2)的数值解?工下面，我们介绍几个常用的边值问题线性多步法，其推导及相关理论可参见Brugnano和Trigiante的专著。三阶广义向后微分公式(GBDF)(12%+1+3%6%_1+yn-2)=h f i,z =2,3,.,A-1,?/3 +6y2 3Vl 2yo)=

27、hfi,|(11W-18W-i+9VN-2 2JN_3)=hfN.五阶广义Adams方法(GAM).一统一1 =4(一19-2+346力_L+456力-74，+i+H+2),=2,.N-2yi-yo=含(19/4+106/3 264/2+646/1+251/0),yN-i UN-2=悬(19加+346/TV-I+456力v2 74册 3 +H/TV-4),、yN VNI=含(251/N+646/-i 264%2 +106/N-3 19/N一4).四阶扩展的第二类梯形方法(ETR2)(盍(%+1+9%2)=亨(力+力i=2,，N -1,击(一 期 +9y2+9 阴17 次)=亨(3一 +/o),

28、点(N3%N2 -9UNI+17g可)=亨(3加 i +于6 六阶方法(TOM)4(ll%+i+2 7yL 2 7 t/i-i 11%2)=(力+i+9力 +9/i-i+力-2),i =2,.,7V 1,tol)for j=1 :N+1Y(n*(j-l)+l:n*j,l)=y(j;F(n*(j l)+l:n*j,:)=feval(rfun,tme sh(j),y(j,:)5);J(n*(j-1)+1:n*j,n*(j l)+l:n*j)=feval(prfun,tme sh(j),y(jendL(1 :N*n,:)=kron(A,I)h*kron(B,I)*J;r(1:N*n)=kron(A,

29、I)*Y+h*kron(B,I)*F;L(n*N+l:n*(N+l),l:n)=fe v al(p g fu n l,y(l,:)?,y(N+lL(n*N+l:n*(N+1),n*N+l:n*(N+l)=feval(pgfun2,y(1,:)5,y(N+1 ,:)5);r(n*N+1:n*(N+1)=fe val(gfun,y(l ,y(N+l xi=Lr;for k=1 :N+1y(k,:)=y(k,:)+xi(n*(k l)+l:n*k)endend end为测试其算法的有效性，我们分别取步长九=1/(5*2，)(，=0,1,2),依次应用上述四种计算格式于例7.2中的边值问题，其数值解的

30、全局误差err(K):=0111a (幻Uj和收敛阶1 r err(h)一P -=l o g 2 777Le rr(/z/2)J列于表7.3中。由该计算结果可见，其算法关于非线性边值问题是非常有效的。hGBDF-3ETR2-4GAM-5TOM-6errhperr Qi)perr(h)perr(h)p1/56.3674e-4-3.1464e-5-1.3510e-5-7.5878e-6-1/109.7062e-52.71372.0233e-63.95894.9753e-74.76314.4113e-87.42631/201.4396e-52.75321.2918e-73.96931.7827e-8

31、4.80269.7162e-105.5047表7.3线性多步法解例7.2的数值结果.7.2 Ritz-Gale rkin方法源于力学中的变分原理有着重要的理论和实践价值，变分原理主要有极小位能原理和虚功原理。在特定的函数空间设置下，微分方程边值问题可转化为与之等价的变分问题，通过求解相应的变分问题，我们可获得原微分方程边值问题的解。其中，基于极小位能原理产生的方法称为Ritz方法，而基于虚功原理产生的方法则称为Gale rkin方法，有时人们也笼统地称这二种方法为Ritz-Gale rkin 方法。7.3.1抽象空间为了方便本节及后面第九章的应用，这里我们首先介绍一些抽象空间概念和常用的积分不

32、等式。设有区域Q u 及巴9Q为其光滑边界，c =为其闭包，则我们可给出如下抽象空间概念及其常用结论：0软。)空间：在。上无穷次可微，在3Q的某邻域上恒为零的全体函数。C7(0)空间：在。上加次连续可微，在3Q的某邻域上恒为零的全体函数。U(0)空间(1 p 00)：在。上夕次绝对可积的全体函数。易知,空间L0)=/:|/|LP(Q)00是一个B a na c h空间。当p =2时，L2(0)空间中的内积和范数可分别定义为此外，在LP(Q)空间中，下列不等式成立：Minkowski 本等式:设/,g L/?(Q),则1 1/+g|i L?(o)-I I/I|LP(Q)+Holde r 不等式：

33、设f e 叫0),g e V(Q),+1=1,贝lj 方 GV(Q),且l l/g|u(o)W -(出I l g l Lg).当夕=q=2时，Holde r不等式即为Cauchy-Schwarz不等式l l/g|i M(c)叫0)。一阶广义导数：设/C L 2(Q),若存在g e L 2(Q),使得g(c)砂(为加=一 j/(劣)袈加，V矽eC软C),JQJ.c，xi其中冗=(的,叼,，喝)丁，则称g是/关于的一阶广义导数，且记为小/(力)。侬阶广义导数:设/6 L2(Q),若存在g 6 L2(Q),使得 g(x)x)d(j=(一1)f(x)dad(y,V 协 e C 1(Q),JoJo其中=

34、(11,比1,，n)T,且Q=(1,Q2,an)e r 同=%,a冲=.端J=I 1则称q为/的|a|阶广义导数，且记为显然广义导数有别于经典导数，前者描述的是全局性质，而后者描述的是局部性质。二者之间有如下关系：若/eC 3(0),则f的广义导数与经典导数相同。例7.3设。=(-1,1),求函数fx=团3 G Q)的一阶广义导数。解V。e Cg(Q),由分部积分可得,IPOP1/(x)d x =/(x)d x +/(x)d xJ-l J-l Jo/0 pl pl=/叭x)dx /叭x)dx=g(x)3(x)d*,J-i Jo J-i其中0-J l，。/11-1,-1 r r 0.因此，根据广

35、义导数的定义可知：d1f=g,但其经典导数不存在。在广义导数的基础上，我们可定义加阶S o b o l e v空 间 间叫。)如下:Hm(Q)=/:VQ:|Q|m,d f e L2(Q),其内积与范数可定义为：(/,g E(加/即匕)近加 n/i im=.J。a po：=mi。M1),q(i)N o(%eQ),f e L2(Q)O对于满足上述边界条件且在。上二次连续可微的任意两函数以力)，N(7),利用分部积分公式可得(Ly,2)=/-五 O)+qy zdx=rbI(pyiz+qyz)dx.a记上式右端为Qg,N),即有Ly,z)=a(y,z).(28)泛函Q,z)将在我们建立边值问题(27)

36、的变分原理时起着举足轻重的作用，特别泛函Q,N)还具有如下重要性质：1)线性性：若 仅。为任意实数，则Q(ayi+眺,z)=aa(gi,z)+(3a(y2,z),Vi/i,y2,z H Q),Q(g,aN i+/?Z2)=a a(g,z i)+/3Q(,Z 2),V g,为，勿 G 卸1。)2)对称性：a(y,z)=a(z,y),Vg,z W 序 3)连续性：存在一个常数M 0,使得|Q(/)|0使得a(y,y)yyt Vg e 畦(。)这里 HI(Q):=y:y e 皿(。)，y(a)=0。显然(Q)C 畦 C 加(Q).连续性的证明：由Cauchy-Schwarz不等式有rbau,v)/(

37、pyrzr+qyz)dxJ aM(J.yzdx+yzd x(=牒 忸(叫,应(叫)故连续性获证。强制性的证明：V ge I(0)有J a因此，由Cauchy-Schwarz不等式得rb、I y(x)dx arba)dx/J a从而rb/yx)2dxJ a这里7=2ay(t)dt dx 7|?/I l i，aa而 端，山。故pb pbQ(y,y)=p(y+qy2dx pQ y2dx yylJ a J a7=Po 7-在泛函Q(g,z)的基础上，我们进一步引入一个新的泛函J(v)=一(4)，ye HN。).该量在物理学中常常表示位能，其使得边值问题(27)的求解问题转化为相应的变分问题。籍此，我们

38、可建立下列变分原理：极小位能原理。定理5.1设函数为 G C2(Q),则以是边值问题(27)的解当且仅当J(y*)=min J(y).(29)证 明Vze HlX。)，我们有Q(g*,z)(/,Z)=/(py*z+qy*z-fz)dxJ a=(py0+/1 一 3pul)+qy*一 f zdx(30)a J a L dx J=p(b)y*(b)z(b)+(Ly*f)zdx.J a若g*是边值问题(27)的解，则有Ly*f =0,n*(b)-o.从而，据(30)得Q(g*,z)(z)=o,V ze H y。).(31)由上式及Q z)的强制性可知：VA G肽0及Vz e HI(Q)0,泛函J)满

39、足Jy*+M =J(y*)+歹。n)J(y*).(32)因此，J(y*)=min J(g)。?/刑(。)反之，若JQ/*)=min J(y),则当记昨明(。)(入)：=J(g*+-z)=J(y*)+Aa(,z)-(/,z)+a(z,z),A e R时，由(30)我们有(0)=P/z(b)+1yLy*-f)zdx=0,V/G HO),(33)且因此/(Ly,-f)zdx=0,/z e Cg(Q).(34)J a应用变分法基本引理于(34)得Ly*一/三 0,/x e Q.(35)将(35)代回(33)得p(b)yz(b)=0,/z e而P(b)Po 0 且取n(2)=x Q时有z(b)O o 因

40、此暝(b)=0o综上所述，以为边值问题(27)的解。该定理要求解外二次连续可微，此时我们称之为边值问题(27)的经典解。而就变分问题(29)而言，其解价不必满足如此强的光滑条件，因此我们称变分问题(29)的解为边值问题(27)的广义解或弱解。以下，我们基于极小位能原理，通过求解变分问题(29)而获得原边值问题(27)的逼近解。设4是回卷(0)的加维子空间(称之力试探函数空 吵：隘是乙 的一组基函数，则/中任一兀素为7)可表mVm=(36)z=l今选取系数匕,使 J(而)取极小。考虑到 7 7 2 mJ(Um)=/(Vm,Vm)(/,淅)=矶外灯)。)，=1 j=l是 4 的二次函数，且Q(2,

41、仍)具对称性，则令-a J(y m)n n 0d c=0,j =1,2,m可得m a(A,0)G =(/,%)，j=1,2,，m.(3 7)2 =1解此关于&的 线 性 方 程 组，且将其解代回(3 6)即得变分问题(29)的逼近解心，这就是R i t z 方法。7.3.3 Gale rkin方3去由 定 理7.1的证明过程，我们可知：若y 02(。)是边值问题(27)的解,则有a(y,n)-(/,z)=0,/z e 卸出。).(38)反之，若V 2(。)是变分方程(38)的解，贝岫(30)得p(b)y(b)z(b)+f Ly f)zdx=0,/z 6J a利用该方程及定理7.1后半部分的证明

42、可推得g是边值问题(27)的解。因此，我们有以下变分原理。定 理7.2设 函 数4 e 2(0),贝切是边值问题(27)的解当且仅当是变分方程(38)的解。由于变分方程(38)的左端在物理学中常常表示虚功，因此我们称 定 理 7.2为虚功原理，而其方程的解也称为边值问题(27)的广义解或弱解。以 下，我 们 基 于 虚 功 原 理，通过求解变分方程(38)而获得原边值问题(27)的逼近解。设%是间玄。)的6 维子空间，6 之是心,的一组基函数，则Um中任一元素而可表示为(36)。今选取系数 Q,使得淅(7)满足变分方程Q(ym,z)=(/,z),/z e Um.(39)m记2=E d用(Vdj

43、 R),将该式及(36)代入(39)得j=im m m dj OQ(p0 0,q(%)0(rr e Q),/G L2(Q)0若记rbQ(g,z)=(pyz+ryz+qyz)dx,J a则我们可建立(40)的虚功原理，并由此导出其Gale rkin方法。但是，由于此时的泛函QG,Z)不再具有对称性及强制性，因此不可能建立(40)的极小位能原理，故而也不可能导出其Ritz方法。例 7.4 取基函数份=s i n(27rrr)(i =1,2,m),用G a l e rk i n方法求解非齐次边值问题(y +=力2,x e (0,1),g(0)=0,g =1.解令如=%u=y y 0,则欲解的非齐次边

44、值问题可化为齐次边值问题(u(x)+(力)=X2 xy x (0,1),(0)=0,M =0.m设该问题的疗汰G a l e rk i n逼近解为明=s，则其满2=1足 R i t z-G a l e rk i n 方 程m。(外均)ci=(/-应日),j=1,2,，mi=l这里Q(A，%)=+(/)%(刈 丘=Jo1(Z 7T)2+1,i=j,0,2平力(x2 x,朽)=/(x2 x)ipj(x)dxJo(1 1.因此,=w 瑞 洲a p L（汨）3 （*产+1,）为奇数,0,2为偶数，且I 7几+1 I8sin(2A:1)TVX昕 :M 侬-1)中(2 k-1)巾+iy这里用表示取整函数。

45、故原边值问题的Gale rkin逼近解为I m+1 I8sin(2k l)7nr而 =力 +)(2 1)中(21闭2+)Ritz-Gale rkin方法导出的逼近解而具有如下收敛性质。定 理7.3设4为试探函数空间，ym e Um是Ritz-Gale rkin方程(39)的解，y e股玄。)是变分方程(38)的解，则存在与y m无关的常数0 0,使得下列误差估计成立：y-ym i 007.3.4有限元法边值问题的Ritz-Gale rkin方 法 存 在 着 若 干 弊 病，首 先，当形成Ritz-Gale rkin方 程 时，我 们 需 要 计 算 大 量 的 积 分，例取“维 试 探 函

46、数 空 间 则 需 要 计 算+3）个（对称情形）到回6+1）个（非对称情形）积分；其次，若Ritz-Gale rkin方程是通过取传统基函数形成的，则该方程组往往是非稀疏的大型病态方程组，其求解十分困难。为克服其困难性，本节介绍有限元方法，该方法与Ritz-Gale rkin方法主要区别在于其选取局部基函数的技巧。目前，有限元方法已被广泛应用于结构力学、固体力学和流体力学等工程领域，其求解的大致步骤如下：S t e p 1.将边值问题转化为变分问题;S t e p 2.将求解域剖分成适当形状的单元;S t e p 3.构造单元形状函数,形成有限元空间;S t e p 4.形成有限元方程,并采

47、用有效方法解之;S t e p 5.分析有限元解。依据上述步骤，我们应用有限元方法于边值问题(2 7),其有限元方程由Gale rkin方法导出。首先，将 求 解 区 间“作网格剖分，网格点为a XQ Xi Xi xm=b,其中小区间Q=称为第。个单元，其长度为居=Xi-Xi-io 此外，我们记最大单元长度/z=max hi，网格点以寸应lim的逼近值为为(i=0,1,m)o取Sobole v空间回卷(0)的一 子 空间Uh=yh：Vh e 11rH(0)n 眶(0),i i m,其 中 为 区 间 0上次数不超过m,的多项式全体，办称为试探函数空间，U/7中的元素仍称为试探函数。以下我们构造

48、线性插值函数构成的试探函数空间，首先我们在每个单元0,上取线性插值y g)=i e 储，(41)由此构成整个求解区间。上的试探函数，这类分段线性函数的全体即构成一个试探函数空间U/,。由(41),试探函数%可在求解区间。上写成统一的形式mzO)=步3,(力)，/e Q,(42)2=1其中 1 XXi,夕?3 =_ N 4?X i X +1,2 =1,2,T儿+1 2.(43)由(37)得有限元方程m a砒 外 伤)%=/fPjdx,J=1,2,m,(44)i=i Ja其中a(,3j)=/(p(pfj H+q(pi(P j)dx,且根据式(43)可知线性方J a程组(44)的系数阵为三对角对称阵

49、/a(夕 1,夕 1)a(2,3 i)Q(2)Q(.2 4 2)Q(夕3,32)K=.a(3 m 2,3m i)a(3 m 1,W m-l)W m-l)a(3 7721,夕m J a(3夕m)/上述有限元方程的矩阵形式为K Y=d,(45)这里y =(阴,V2,m)T,d=(di,d2,dm)T,4=/f(pidx,J aK、！称R X 为m总 A刚A e度 H矩R -阵，三。、利U&7用J追赶法解此方程组，我们即可获得边值例7.5用有限元法解边值问题(y+?/=(1 +7r2)sin(7r1),x e(0,1),g(0)=0,g=0.解 将 0,1 作m等分，记 九=2，Xi=Xi-i+ho

50、既然边界值为已知，故我们只需取巾-1维试探函数空间U/z，其基函数可取为(X Xi-I7,*1 X X i,h仍(7)=11 _ xX l xx x 1，/=1,2,，加 L7 )：儿 _ 4 2+1)、0,oth erwise因此，有限元解可表示为m 1yh(x)=f y 捋 3,2=1其中u=%a4通过追赶法解有限元方程a（3 2,初）Q（6 2）。（夕2,3 2）。（夕3 2）aSm-3）37n2）。（夕7n2,Q（份 九 2,3772 2)a(夕m,1,3m 2)夕m 1)。(夕m 1,3m 1,1、V2ym-i/ddmi获得，这里a ,%）=/（d *+6（Pj）dx,di=其有限元