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

偏微分方程数值解上机实验报告（一）实验一上机题目：用线性元求解下列边值问题的数值解：-y''+π24y=π22sinπ2x0<x<1y(0)=0y'1=0实验程序：functionS=bzx=fzero(@zfun1);[ty]=ode45(@odefun[01][0x]);S.t=t;S.y=y;plot(ty)xlabel('x:´从0一直到1')ylabel('y')title('线性元求解边值问题的数值解')functiondy=odefun(xy)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<x<1区间上随着x的不断变化xy之间关系如下图所示：（二）实验二上机题目：求解Helmholtz方程的边值问题：于于其中k=15101520实验程序：（采用有限元方法这里对[01]*[01]采用n*n的划分n为偶数）n=10;a=zeros(n);f=zeros(n);b=zeros(1n);U=zeros(n1);u=zeros(n1);fori=2:na(i-1i-1)=pi^2/(12*n)+n;a(i-1i)=pi^2/(24*n)-n;a(ii-1)=pi^2/(24*n)-n;forj=1:nifj==i-1a(ij)=a(ii-1);elseifj==ia(i-1j-1)=2*a(i-1i-1);elseifj==i+1a(ij)=a(ii+1);elsea(ij)=0;endendendendenda(nn)=pi^2/(12*n)+n;fori=2:nf(i-1i)=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);endfori=1:nf(ii)=-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-1j)+f(ij)fori=1:(n-1)b(i)=f(ii)+f(ii+1);endb(n)=f(nn);tic;n=20;can=20;s=zeros(n^210);h=1/n;st=1/(2*n^2);A=zeros((n+1)^2(n+1)^2);symsxy;fork=1:1:2*n^2s(k1)=k;q=fix(k/(2*n));r=mod(k(2*n));if(r~=0)r=r;elser=2*n;q=q-1;endif(r<=n)s(k2)=q*(n+1)+r;s(k3)=q*(n+1)+r+1;s(k4)=(q+1)*(n+1)+r+1;s(k5)=(r-1)*h;s(k6)=q*h;s(k7)=r*h;s(k8)=q*h;s(k9)=r*h;s(k10)=(q+1)*h;elses(k2)=q*(n+1)+r-n;s(k3)=(q+1)*(n+1)+r-n+1;s(k4)=(q+1)*(n+1)+r-n;s(k5)=(r-n-1)*h;s(k6)=q*h;s(k7)=(r-n)*h;s(k8)=(q+1)*h;s(k9)=(r-