椭圆型方程的有限元法

两点边值问题有限元法(必做)从Galerkin原理出发用线性元解两点边值问题-u”+u=x2,0<x<1U(0)=U(l)=0精确解：u(x)=―1—[(2一3e)ex一(2e一3)e1-x]+x2+2。e2-11.1变分形式从Galerkin原理出发推导出两点边值问题的变分形式，将积分区间等分为N份’则步长h=1i=口"’记为h。写出有限元方程及系数矩阵元素。基于虚功原理，求变分形式a(u,v)=(f,v)。hx=0,x,…,x=1。取值为01n1<x=0,x,…,x=1。取值为01n1<i<N,h=—。取u=v,

NTOC\o"1-5"\h\zhhhu(0)=u=0,u，…,u=u(1)=0。其中x=x+ih，01ni0J1(-u''+u)udx=J1x2udx，推得A[(u')2+u2]dx=J1ux2dx。相应的双线性变分形式0000a(Q,Q)=J1[Q'Q+QQ]dx，则有限元方程a(Q(x),Q(x))u=(f(x),Q(x))，ij0ijijijiji=1a仰,甲)=J1[-h-1p+hq(1-8)8]d8；j-1j0jja(Q,Q)=J1[h-1p+hq82]d8+J1[h-1q(1-8)2]d8；jj0jj0j+1a(Q+Q)=Jxj+1[—h-1p+hq8(1—8)]d8；j+1jxj+1j+1j这里j=2,…,n-1。第一行只有两个非零元素：a(Q,Q)，a(Q,Q)。第n行1112n-1nnJ*1dx=hJ1f(x+h8)ds+hJ1f(x0jj0j-1n-1nnJ*1dx=hJ1f(x+h8)ds+hJ1f(x0jj0j-1jj+10jj+1nn-1nn0nn

TOC\o"1-5"\h\za(p,p)a(p,p)…01112方程的系数矩阵为:a(p,p)a(p,p)•…0方程的系数矩阵为:212.2.::••・a(p,p)n-1n00…a(p,p)nn1.2利用MATLAB求解问题的过程依次取N二2n,n二2,3,4,5,6,7,8.用MATLAB求解并图形比较数值解与精确解,用表格列出不同剖分时的L2误差。N=4：Figure1口FileEditViewInsertToolsDesktop\Vindaw旦mlp■ZI阴□口N=8:N=16:

N=32：N二64：

N=128：N=256：误差比较：N48163264128256Err0.0002460.0000860.0000300.0000110.0000040.0000010.0000001.3方法总结及分析在利用Galerkin原理出发用线性元解两点边值问题，利用MATLAB作图可以发现解析解与精确解非常逼近，但从误差上可以看出，剖分结点越多，误差越小，逼近程度越好。附件程序function[U,precise_value,err]=G(N)h=1/N;p=1;q=1;X=0:h:1;A=zeros(N-1);fori=2:N-1f3=@(ks)-p./h+h.*q.*ks.*(1-ks);f2=@(ks)p./h+h.*q.*(ks.入2)+p./h+h.*q.*((1—ks).入2);f1=@(ks)-p./h+h.*q.*ks.*(1-ks);A(i-1,i)=quadl(f1,0,1);A(i,i)=quadl(f2,0,1);A(i,i-1)=quadl(f3,0,1);endA(1,1)=quadl(f2,0,1);f=zeros(N-1,1);fori=2:Nf11=@(ks)(X(i—1)+h.*ks).入2.*ks+(X(i)+h.*ks).入2.*(1—ks);f(i-1)=h.*quadl(f11,0,1);endU=A\f;dx=X;precise_value=((exp(2)-1)入(-1)).*((2-3*exp(1))*exp(dx)-(2*exp(1)-3)*exp(1-dx))+dx.入2+2plot(X,[0;U;0],'b--',X,precise_value,'r:+');legend('数值解','精确解');err=norm