偏微分方程数值解实验报告.docx
- 文档编号:8789127
- 上传时间:2023-02-01
- 格式:DOCX
- 页数:11
- 大小:110.33KB
偏微分方程数值解实验报告.docx
《偏微分方程数值解实验报告.docx》由会员分享,可在线阅读,更多相关《偏微分方程数值解实验报告.docx(11页珍藏版)》请在冰豆网上搜索。
偏微分方程数值解实验报告
偏微分方程数值解
上
机
实
验
报
告
(一)实验一
一、上机题目:
用线性元求解下列边值问题的数值解:
-
sin
,0 y(0)=0, 二、实验程序: functionS=bz x=fzero(@zfun,1); [ty]=ode45(@odefun,[01],[0x]); S.t=t; S.y=y; plot(t,y) xlabel('x: ´从0一直到1') ylabel('y') title('线性元求解边值问题的数值解') functiondy=odefun(x,y) dy=[00]'; dy (1)=y (2); dy (2)=(pi^2)/4*y (1)-((pi^2)/2)*sin(x*pi/2); functionz=zfun(x); [ty]=ode45(@odefun,[01],[0x]); z=y(end)-0; 三、实验结果: 1.以步长h=0.05进行逐步运算,运行上面matlab程序结果如下: 2.在0 (二)实验二 四、上机题目: 求解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); fori=2: n a(i-1,i-1)=pi^2/(12*n)+n; a(i-1,i)=pi^2/(24*n)-n; a(i,i-1)=pi^2/(24*n)-n; forj=1: n ifj==i-1 a(i,j)=a(i,i-1); elseifj==i a(i-1,j-1)=2*a(i-1,i-1); elseifj==i+1 a(i,j)=a(i,i+1); else a(i,j)=0; end end end end end a(n,n)=pi^2/(12*n)+n; fori=2: n f(i-1,i)=4/pi*cos((i-1)*pi/2/n)-8*n/(pi^2)*sin(i*pi/2/n)+8*n/(pi^2)*sin((i-1)*pi/2/n); end fori=1: n f(i,i)=-4/pi*cos(i*pi/2/n)+8*n/(pi^2)*sin(i*pi/2/n)-8*n/(pi^2)*sin((i-1)*pi/2/n); end %b(j)=f(i-1,j)+f(i,j) fori=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(n^2,10); h=1/n; st=1/(2*n^2); A=zeros((n+1)^2,(n+1)^2); symsxy; fork=1: 1: 2*n^2 s(k,1)=k; q=fix(k/(2*n)); r=mod(k,(2*n)); if(r~=0) r=r; elser=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); fork=1: 1: 2*n^2 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); fori=1: 1: 3 forj=i: 3 d(i,j)=int(int(((((diff(L(i),x))*(diff(L(j),x)))+((diff(L(i),y))*(diff(L(j),y))))-((can^2)*L(i)*L(j))),x,0,1),y,0,1); d(j,i)=d(i,j); end end fori=1: 1: 3 forj=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 fori=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,n^2); j=n^2; fori=(n^2+n): -1: 1 if((mod(i,(n+1)))~=1) M(: j)=A(: i); j=j-1; elsecontinue end end preanswer=M\B; answer=zeros((n+1)^2,1); j=1; fori=1: 1: (n^2+n) if((mod(i,(n+1)))~=1) answer(i)=preanswer(j); j=j+1; elseanswer(i)=0; end end Z=zeros((n+1),(n+1)); fori=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 fori=1: n u(i,1)=4/(pi^2)*sin(pi*i/n/2); end u e=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的编程语法。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 实验 报告