偏微分方程数值解实验报告1、用有限元方法求下列边值问题的数值解：''()112x -y +y =2sin ,0<x <,42y(0)=0,y'()=0.ππ⎧⎪⎨⎪⎩其中其精确解为 24x y =sin()2ππ，取h=0.1要求：（1）将精确解与用有限元得到的数值解画在同一图中（2）1hH u u -、2h L u u -、01h x max u -u ≤≤2、用线性元求解下列问题的数值解：x x u(x,-1)=u(x,1)=0,-1<x <1,u (-1,y)=1,u u =0=-2,-1<x,,-1<y <1.y <1,⎧∆⎪⎨⎪⎩ 精确到小数点后第六位，并画出解曲面。

3、用Crank-Nicolson 差分法求解Burger 方程0,(0,1),(0,5),,[01]t (1).t x xx v +vv =v ,x t v(x,0)=sin2x x ;v =v ,t =0ννπ>∈∈⎧⎨∈(0,)⎩，其中取1ν= 要求画出解曲面。

迭代格式如下：1221212111111111122142212n nn n n n j j j j j j n n n n n n j j j j j j V V V V V V h h V V V V V V h h τ++++++++++-+-⎧⎫-()-()()-()⎪⎪++⎨⎬⎪⎪⎩⎭⎧⎫-+-+⎪⎪=+⎨⎬⎪⎪⎩⎭1、%Ritz Galerkin方法求解方程function u1=Ritz(x)%定义步长h=1/100;x=0:h:1;n=1/h;a=zeros(n-1,1);b=zeros(n,1);c=zeros(n-1,1);d=zeros(n,1);%求解Ritz方法中内点系数矩阵for i=1:1:n-1b(i)=(1/h+h*pi*pi/12)*2;d(i)=h*pi*pi/2*sin(pi/2*(x(i)+h))/2+h*pi*pi/2*sin(pi/2*x(i+1))/2; end%右侧导数条件边界点的计算b(n)=(1/h+h*pi*pi/12);d(n)=h*pi*pi/2*sin(pi/2*(x(i)+h))/2;for i=1:1:n-1a(i)=-1/h+h*pi*pi/24;c(i)=-1/h+h*pi*pi/24;end%调用追赶法u=yy(a,b,c,d)%得到数值解向量u1=[0,u]%对分段区间做图plot(x,u1)%得到解析解y1=sin(pi/2*x);hold onplot(x,y1,'o')legend('数值解','解析解')function x=yy(a,b,c,d)n=length(b);q=zeros(n,1);p=zeros(n,1);q(1)=b(1);p(1)=d(1);for i=2:1:nq(i)=b(i)-a(i-1)*c(i-1)/q(i-1);p(i)=d(i)-p(i-1)*c(i-1)/q(i-1);endx(n)=p(n)/q(n);for j=n-1:-1:1x(j)=(p(j)-a(j)*x(j+1))/q(j);endxx =Columns 1 through 110.0157 0.0314 0.0471 0.0628 0.0785 0.0941 0.1097 0.1253 0.1409 0.1564 0.1719Columns 12 through 220.1874 0.2028 0.2181 0.2335 0.2487 0.2639 0.2790 0.2940 0.3090 0.3239 0.3387Columns 23 through 330.3535 0.3681 0.3827 0.3972 0.4115 0.4258 0.4400 0.4540 0.4679 0.4818 0.4955Columns 34 through 440.5091 0.5225 0.5358 0.5490 0.5621 0.5750 0.5878 0.6004 0.6129 0.6253 0.6374Columns 45 through 550.6495 0.6613 0.6730 0.6846 0.6959 0.7071 0.7181 0.7290 0.7397 0.7501 0.7604Columns 56 through 660.7705 0.7805 0.7902 0.7997 0.8090 0.8182 0.8271 0.8358 0.8444 0.8527 0.8608Columns 67 through 770.8687 0.8763 0.8838 0.8910 0.8981 0.9049 0.9114 0.9178 0.9239 0.9298 0.9355Columns 78 through 880.9409 0.9461 0.9511 0.9558 0.9603 0.9646 0.9686 0.9724 0.9759 0.9793 0.9823Columns 89 through 990.9851 0.9877 0.9901 0.9921 0.9940 0.9956 0.9969 0.9981 0.9989 0.9995 0.9999Column 1001.0000u =Columns 1 through 110.0157 0.0314 0.0471 0.0628 0.0785 0.0941 0.1097 0.1253 0.1409 0.1564 0.1719Columns 12 through 220.1874 0.2028 0.2181 0.2335 0.2487 0.2639 0.2790 0.2940 0.3090 0.3239 0.3387Columns 23 through 330.3535 0.3681 0.3827 0.3972 0.4115 0.4258 0.4400 0.4540 0.4679 0.4818 0.4955Columns 34 through 440.5091 0.5225 0.5358 0.5490 0.5621 0.5750 0.5878 0.6004 0.6129 0.6253 0.6374Columns 45 through 550.6495 0.6613 0.6730 0.6846 0.6959 0.7071 0.7181 0.7290 0.7397 0.7501 0.7604Columns 56 through 660.7705 0.7805 0.7902 0.7997 0.8090 0.8182 0.8271 0.8358 0.8444 0.8527 0.8608Columns 67 through 770.8687 0.8763 0.8838 0.8910 0.8981 0.9049 0.9114 0.9178 0.9239 0.9298 0.9355Columns 78 through 880.9409 0.9461 0.9511 0.9558 0.9603 0.9646 0.9686 0.9724 0.9759 0.9793 0.9823Columns 89 through 990.9851 0.9877 0.9901 0.9921 0.9940 0.9956 0.9969 0.9981 0.9989 0.9995 0.9999Column 1001.0000u1 =Columns 1 through 110 0.0157 0.0314 0.0471 0.0628 0.0785 0.0941 0.1097 0.1253 0.1409 0.1564Columns 12 through 220.1719 0.1874 0.2028 0.2181 0.2335 0.2487 0.2639 0.2790 0.2940 0.3090 0.3239Columns 23 through 330.3387 0.3535 0.3681 0.3827 0.3972 0.4115 0.4258 0.4400 0.4540 0.4679 0.4818Columns 34 through 440.4955 0.5091 0.5225 0.5358 0.5490 0.5621 0.5750 0.5878 0.6004 0.6129 0.6253Columns 45 through 550.6374 0.6495 0.6613 0.6730 0.6846 0.6959 0.7071 0.7181 0.7290 0.7397 0.7501Columns 56 through 660.7604 0.7705 0.7805 0.7902 0.7997 0.8090 0.8182 0.8271 0.8358 0.8444 0.8527Columns 67 through 770.8608 0.8687 0.8763 0.8838 0.8910 0.8981 0.9049 0.9114 0.9178 0.9239 0.9298Columns 78 through 880.9355 0.9409 0.9461 0.9511 0.9558 0.9603 0.9646 0.9686 0.9724 0.9759 0.9793Columns 89 through 990.9823 0.9851 0.9877 0.9901 0.9921 0.9940 0.99560.9969 0.9981 0.9989 0.9995Columns 100 through 1010.9999 1.0000ans =Columns 1 through 100 0.0157 0.0314 0.0471 0.0628 0.0785 0.0941 0.1097 0.1253 0.1409Columns 11 through 200.1564 0.1719 0.1874 0.2028 0.2181 0.2335 0.2487 0.2639 0.2790 0.2940Columns 21 through 300.3090 0.3239 0.3387 0.3535 0.3681 0.3827 0.3972 0.4115 0.4258 0.4400Columns 31 through 400.4540 0.4679 0.4818 0.4955 0.5091 0.5225 0.5358 0.5490 0.5621 0.5750Columns 41 through 500.5878 0.6004 0.6129 0.6253 0.6374 0.6495 0.6613 0.6730 0.6846 0.6959Columns 51 through 600.7071 0.7181 0.7290 0.7397 0.7501 0.7604 0.7705 0.7805 0.7902 0.7997Columns 61 through 700.8090 0.8182 0.8271 0.8358 0.8444 0.8527 0.8608 0.8687 0.8763 0.8838Columns 71 through 800.8910 0.8981 0.9049 0.9114 0.9178 0.9239 0.9298 0.9355 0.9409 0.9461Columns 81 through 900.9511 0.9558 0.9603 0.9646 0.9686 0.9724 0.9759 0.9793 0.9823 0.9851Columns 91 through 1000.9877 0.9901 0.9921 0.9940 0.9956 0.9969 0.9981 0.9989 0.9995 0.9999Column 1011.00002、function [ u ] = Q_2( P )format longif nargin<1P=16;endf=2;beta=-1;x=linspace(-1,1,P+1);h=2/P;%%%Q=P*P-1;Qi=(P-1)*(P-1);pos=zeros(Q,2);i=2;j=2;for n=1:Qipos(n,1)=x(i);pos(n,2)=x(j);j=j+1;if mod(n,P-1)==0i=i+1;j=2;endendfor j=1:P-1pos(j+Qi,1)=-1;pos(j+Qi,2)=x(j+1);pos(j+Qi+P-1,1)=1;pos(j+Qi+P-1,2)=x(j+1);end%%b=zeros(Q,1);for n=(Qi+1):(Qi+P-1)b(n)=beta*h;endfv=f*h*h/6;for n=1:Qib(n)=b(n)+6*fv;endfor n=1:P-1b(Qi+n)=b(Qi+n)+fv*3;b(Qi+n+P-1)=b(Qi+n+P-1)+fv*3;end%%K=zeros(Q);for n=1:QiK(n,n)=4;endfor n=Qi+1:Qi+P-1K(n,n)=2;K(n+P-1,n+P-1)=2;endfor i=1:Qifor j=i+1:Qif abs(pos(i,1)-pos(j,1))+abs(pos(i,2)-pos(j,2))==hK(i,j)=-1;end;endendfor i=Qi+1:Qi+P-1for j=i+1:Qi+P-1if abs(pos(i,1)-pos(j,1))+abs(pos(i,2)-pos(j,2))==hK(i,j)=-0.5;end;if abs(pos(i+P-1,1)-pos(j+P-1,1))+abs(pos(i+P-1,2)-pos(j+P-1,2))==h K(i+P-1,j+P-1)=-0.5;endendendfor i=1:Qfor j=i+1:QK(j,i)=K(i,j);endendu=linsolve(K,b);[X,Y]=meshgrid(x,x);U=zeros(P+1);U(2:P,1)=u(Qi+1:Qi+P-1);U(2:P,P+1)=u(Qi+P:Q);U(2:P,2:P)=reshape(u(1:Qi),P-1,P-1);surf(X,Y,U);xlabel('x');ylabel('y');zlabel('u');str=strcat(['N =',num2str(P)],[', 步长h =',num2str(h)]); title(str);hold offend3、a=0;b=1;t1=0;t2=5;h=0.01;tao=0.01;n1=(t2-t1)/tao;n2=(b-a)/h;eps=10^-8;V=zeros(n1+1,n2+1);M=zeros(n2-1,n2-1);m=zeros(n2-1,1);x=zeros(n2+1,1);y=zeros(n2+1,1);y1=zeros(n2+1,1);temp=zeros(n2+1,1);for i=2:n2r=(i-1)*h;V(1,i)=sin(2*pi*r);x(i)=V(1,i);enda=1/tao+1/h^2;y=x;for i=2:n1+1while norm((y-temp),2)>eps %应改为向量2-范数%构造Mfor j=1:n2-1if j==1M(1,1)=a-y(2);M(1,2)=y(3)/(2*h)-1/(2*h^2);endif j==n2-1M(n2-1,n2-2)=-1/(2*h^2);M(n2-1,n2-1)=a-y(n2)/(2*h);endif j~=1 && j~=n2-1M(j,j-1)=-1/(2*h^2);M(j,j)=a-y(j+1);M(j,j+1)=y(j+2)/(2*h)-1/(2*h^2);endend%构造mfor j=1:n2-1m(j)=(y(j+1)-x(j+1))/tao+(x(j+2)^2-x(j+1)^2+y(j+2)^2-y(j+1)^2)/(4*h)-(x(j+2)-2*x(j+1)+x(j)+y(j+2)-2*y(j+1)+y(j))/(2*h^2);endtemp=y;y(2:n2)=y(2:n2)-inv(M)*m;%y=y1;endV(i,:)=y';x=y;y=zeros(n2+1,1);endxx = 0:h:1;yy = 0:tao:5;surf(xx,yy,V)。