偏微分方程数值解实验报告.doc

偏微分方程数值解上机实验报告 （一）实验一一、 上机题目：用线性元求解下列边值问题的数值解：-y+24y=22sin2x，0x1y(0)=0，y1=0二、 实验程序：function S=bz x=fzero(zfun,1); t y=ode45(odefun,0 1,0 x); S.t=t;S.y=y;plot(t,y)xlabel(x:从0一直到1)ylabel(y) title(线性元求解边值问题的数值解) function dy=odefun(x,y) dy=0 0;dy(1)=y(2);dy(2)=(pi2)/4*y(1)-(pi2)/2)*sin(x*pi/2);function z=zfun(x); t y=ode45(odefun,0 1,0 x);z=y(end)-0;三、 实验结果：1以步长h=0.05进行逐步运算，运行上面matlab程序结果如下： 2在0x1区间上，随着x的不断变化，x，y之间关系如下图所示：（二）实验二四、 上机题目：求解Helmholtz方程的边值问题：，于，于其中k=1,5,10,15,20五、 实验程序：（采用有限元方法，这里对0,1*0,1采用n*n的划分，n为偶数）n=10; a=zeros(n);f=zeros(n);b=zeros(1,n);U=zeros(n,1);u=zeros(n,1); for i=2:n a(i-1,i-1)=pi2/(12*n)+n; a(i-1,i)= pi2/(24*n)-n; a(i,i-1)= pi2/(24*n)-n; for j=1:n if j=i-1 a(i,j)=a(i,i-1); else if j=i a(i-1,j-1)=2*a(i-1,i-1); else if j=i+1 a(i,j)=a(i,i+1); else a(i,j)=0; end end endendenda(n,n)=pi2/(12*n)+n; for i=2:n f(i-1,i)=4/pi*cos(i-1)*pi/2/n)-8*n/(pi2)*sin(i*pi/2/n)+8*n/(pi2)*sin(i-1)*pi/2/n); end for i=1:n f(i,i)=-4/pi*cos(i*pi/2/n)+8*n/(pi2)*sin(i*pi/2/n)-8*n/(pi2)*sin(i-1)*pi/2/n); end %b(j)=f(i-1,j)+f(i,j)for i=1:(n-1) b(i)=f(i,i)+f(i,i+1); end b(n)=f(n,n); tic; n=20; can=20; s=zeros(n2,10); h=1/n; st=1/(2*n2); A=zeros(n+1)2,(n+1)2); syms x y; for k=1:1:2*n2 s(k,1)=k; q=fix(k/(2*n); r=mod(k,(2*n); if (r=0) r=r; else r=2*n;q=q-1; end if (r=n) s(k,2)=q*(n+1)+r; s(k,3)=q*(n+1)+r+1; s(k,4)=(q+1)*(n+1)+r+1; s(k,5)=(r-1)*h; s(k,6)=q*h; s(k,7)=r*h; s(k,8)=q*h; s(k,9)=r*h; s(k,10)=(q+1)*h; else s(k,2)=q*(n+1)+r-n; s(k,3)=(q+1)*(n+1)+r-n+1; s(k,4)=(q+1)*(n+1)+r-n; s(k,5)=(r-n-1)*h; s(k,6)=q*h; s(k,7)=(r-n)*h; s(k,8)=(q+1)*h; s(k,9)=(r-n-1)*h; s(k,10)=(q+1)*h; end end d=zeros(3,3); B=zeros(n+1)2,1); b=zeros(3,1); for k=1:1:2*n2 L(1)=(1/(2*st)*(s(k,7)*s(k,10)-s(k,9)*s(k,8)+(s(k,8)-s(k,10)*x+(s(k,9)-s(k,7)*y); L(2)=(1/(2*st)*(s(k,9)*s(k,6)-s(k,5)*s(k,10)+(s(k,10)-s(k,6)*x+(s(k,5)-s(k,9)*y); L(3)=(1/(2*st)*(s(k,5)*s(k,8)-s(k,7)*s(k,6)+(s(k,6)-s(k,8)*x+(s(k,7)-s(k,5)*y); for i=1:1:3 for j=i:3 d(i,j)=int(int(diff(L(i),x)*(diff(L(j),x)+(diff(L(i),y)*(diff(L(j),y)-(can2)*L(i)*L(j),x,0,1),y,0,1); d(j,i)=d(i,j); end end for i=1:1:3 for j=1:1:3 A(s(k,(i+1),s(k,(j+1)=A(s(k,(i+1),s(k,(j+1)+d(i,j); end end for i=1:1:3 b(i)=int(int(L(i),x,0,1),y,0,1); B(s(k,(i+1),1)=B(s(k,(i+1),1)+b(i); end end M=zeros(n+1)2,n2); j=n2; for i=(n2+n):-1:1 if (mod(i,(n+1)=1) M(:,j)=A(:,i); j=j-1; else continue end end preanswer=MB; answer=zeros(n+1)2,1); j=1; for i=1:1:(n2+n) if (mod(i,(n+1)=1) answer(i)=preanswer(j); j=j+1; else answer(i)=0; end end Z=zeros(n+1),(n+1); for i=1:1:(n+1)2 s=fix(i/(n+1)+1; r=mod(i,(n+1); if(r=0) r=n+1; s=s-1; else end Z(r,s)=answer(i); end X,Y=meshgrid(1:-h:0,0:h:1); surf(X,Y,Z); toc; t=toc; K=a ;B=b; U=inv(K)*B for i=1:n u(i,1)=4/(pi2)*sin(pi*i/n/2);end ue=U-u六、 实验结果：程序中的变量can为问题中的k，为了方便比较，采用了画图的方式。在n=10时，分别考察can=k=1,5,10,15,20时的情况。1.can=k=1时2.can=k=5时3.can=k=10时4.can=k=15时5.can=k=20时六、实验总结 本次实验第二题很难，通过查找资料和请教同学，学到了很多。通过本次实验，我也进一步了解了线性元求边值问题的数值解的思想和方法，同时了解了matlab工具包中相关类库对求边值问题数值解的支持，更加熟悉了matlab的编程语法。