偏微分方程数值实验报告六实验题目：用Ritz-Galerkin 方法求边值问题2(0)0,(1)1u u x u u ''-+===⎧⎨⎩01x << 的第n 次近似()n u x ，基函数为(),i 1,2,.()..x n i sin i x πφ==并用表格列出0.25，0.5,0.75三点处的真解和1,2,3,4n =时的数值解。

实验算法：将上述边值问题转化为基于虚功方程的变分问题为： 求10()u H I ∈，满足10(,)(,(v H ,))a u v x x I ∀∈=-其中1100(u''v uv)(,)(,)dx (u'v'uv)dx a u v Lu v +=-+==⎰⎰ 记(x)sin()w x π=，引入10()H I 的n 维近似子空间12{,},,n n U φφφ=⋯(),i 1,2,,i sin i x n πφ==⋯利用Ritz-Galerkin 公式：1()(f,),j 1,,,2,n j i jj a n φφφ==⋯=∑，则问题关于n U 下的近似变分问题解1()()n i i n i c u x x φ==∑中的系数12(,,,)T nn c c c c =⋯∈满足12()(,,)n j i j j a x φφφ==∑程序代码：%主函数pde=modeldata();I=[0,1];%积分精度quadorder=10;n=[1,2,3,4];for i=1:length(n)uh{i}=Galerkin(pde,I,n(i),quadorder); endshowsolution(uh{1},'-bo');hold onshowsolution(uh{2},'-rx');hold onshowsolution(uh{3},'-.ko');hold onshowsolution(uh{4},'--k*');hold offtitle('the solution of the un');xlabel('x');ylabel('u');legend('n=1','n=2','n=3','n=4');%计算这三点的数值解x=[1/4,1/2,3/4];un=zeros(length(n),3);for i=1:length(n)[v, ]=basis(x,n(i));un(i,:)=(v'*uh{i})';endformat shortedisp('un');disp(un);%存储数据function pde=modeldata()pde=struct('f',@f);function z=f(x)z=x.*x;endend%Galerkin方法function uh=Galerkin(pde,I,n,quadorder) h=I(2)-I(1);[lambda,weight]=quadpts1d(I,quadorder);nQuad=length(weight);A=zeros(n,n);b=zeros(n,1);for q=1:nQuadgx=lambda(q);w=weight(q);x=[0.25;0.5;0.75];[phi,gradPhi]=basis(gx,n);A=A+(-gradPhi*gradPhi'+phi*phi')*w; b=b+pde.f(gx)*phi*w;endA=h*A;b=h*b;uh=A\b;end%构造基函数function [phi,gradPhi]=basis(x,n)m=length(x);w=sin(pi*x);v=ones(n,m);v(2:end,:)=bsxfun(@times,v(2:end,:),x);v=cumprod(v,1);phi=bsxfun(@times,v,w);gw=pi*cos(pi*x);gv=[zeros(1,m);v(1:end-1,:)];gv(3:end,:)=bsxfun(@times,(2:n-1)',gv(3:end,:)); gradPhi=bsxfun(@times,v,gw)+bsxfun(@times,gv,w); end%数值解的图像function showsolution(u,varargin)x=0:0.01:1;n=length(u);[v, ]=basis(x,n);y=v'*u;plot(x,y,varargin{:});% varargin: an input variable in a functionend%系数矩阵Afunction [lambda,weight] = quadpts1d(I,quadorder)numPts = ceil((quadorder+1)/2);if numPts> 10numPts = 10;endswitch numPtscase 1A = [0 2.0000000000000000000000000];case 2A = [0.57735026918962576450914881.0000000000000000000000000-0.57735026918962576450914881.0000000000000000000000000];case 3A = [ 00.88888888888888888888888890.77459666924148337703585310.5555555555555555555555556-0.7745966692414833770358531 0.5555555555555555555555556];case 4A = [0.33998104358485626480266580.65214515486254614262693610.8611363115940525752239465 0.3478548451374538573730639-0.3399810435848562648026658 0.6521451548625461426269361-0.8611363115940525752239465 0.3478548451374538573730639];case 5A = [ 00.56888888888888888888888890.5384693101056830910363144 0.47862867049936646804129150.9061798459386639927976269 0.2369268850561890875142640-0.5384693101056830910363144 0.4786286704993664680412915-0.9061798459386639927976269 0.2369268850561890875142640];case 6A = [0.23861918608319690863050170.46791393457269104738987030.6612093864662645136613996 0.36076157304813860756983350.9324695142031520278123016 0.1713244923791703450402961-0.2386191860831969086305017 0.4679139345726910473898703-0.6612093864662645136613996 0.3607615730481386075698335-0.9324695142031520278123016 0.1713244923791703450402961];case 7A = [00.41795918367346938775510200.4058451513773971669066064 0.38183005050511894495036980.7415311855993944398638648 0.27970539148927666790146780.9491079123427585245261897 0.1294849661688696932706114-0.4058451513773971669066064 0.3818300505051189449503698-0.7415311855993944398638648 0.2797053914892766679014678-0.9491079123427585245261897 0.1294849661688696932706114];case 8A = [0.18343464249564980493947610.36268378337836198296515040.5255324099163289858177390 0.31370664587788728733796220.7966664774136267395915539 0.22238103445337447054435600.9602898564975362316835609 0.1012285362903762591525314-0.1834346424956498049394761 0.3626837833783619829651504-0.5255324099163289858177390 0.3137066458778872873379622-0.7966664774136267395915539 0.2223810344533744705443560-0.9602898564975362316835609 0.1012285362903762591525314];case 9A = [00.33023935500125976316452510.3242534234038089290385380 0.31234707704000284006863040.6133714327005903973087020 0.26061069640293546231874290.8360311073266357942994298 0.18064816069485740405847200.9681602395076260898355762 0.0812743883615744119718922-0.3242534234038089290385380 0.3123470770400028400686304-0.6133714327005903973087020 0.2606106964029354623187429-0.8360311073266357942994298 0.1806481606948574040584720-0.9681602395076260898355762 0.0812743883615744119718922];case 10A = [0.14887433898163121088482600.29552422471475287017389300.4333953941292471907992659 0.26926671930999635509122690.6794095682990244062343274 0.21908636251598204399553490.8650633666889845107320967 0.14945134915058059314577630.9739065285171717200779640 0.0666713443086881375935688-0.1488743389816312108848260 0.2955242247147528701738930-0.4333953941292471907992659 0.2692667193099963550912269-0.6794095682990244062343274 0.2190863625159820439955349-0.86506336668898451073209670.1494513491505805931457763-0.97390652851717172007796400.0666713443086881375935688];endlambda = (I(2)+I(1)+(I(2)-I(1))*A(:,1))/2; weight = A(:,2)*(I(2)-I(1))/2;数值解：x=x=3/4x=1/21/4u-3.0184e-02 -4.2686e-02 -3.0184e-02 1u-2.1919e-02 -4.2686e-02 -3.8448e-02 2u-2.3232e-02 -4.0652e-02 -3.9761e-02 3u-2.3472e-02 -4.0652e-02 -3.9521e-02 4图像：00.10.20.30.40.50.60.70.80.91the solution of the unx u。