橢圓型方程的有限元法

兩點邊值問題有限元法(必做)從Galerkin原理出發(fā)用線性元解兩點邊值問題-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原理出發(fā)推導(dǎo)出兩點邊值問題的變分形式，將積分區(qū)間等分為N份’則步長h=1i=口"’記為h。寫出有限元方程及系數(shù)矩陣元素?；谔摴υ?，求變分形式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。相應(yīng)的雙線性變分形式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方程的系數(shù)矩陣為:a(p,p)a(p,p)?…0方程的系數(shù)矩陣為:212.2.::???a(p,p)n-1n00…a(p,p)nn1.2利用MATLAB求解問題的過程依次取N二2n,n二2,3,4,5,6,7,8.用MATLAB求解并圖形比較數(shù)值解與精確解,用表格列出不同剖分時的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方法總結(jié)及分析在利用Galerkin原理出發(fā)用線性元解兩點邊值問題，利用MATLAB作圖可以發(fā)現(xiàn)解析解與精確解非常逼近，但從誤差上可以看出，剖分結(jié)點越多，誤差越小，逼近程度越好。附件程序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('數(shù)值解','精確解');err=norm