(完整word版)偏微分方程数值解实验报告
偏徽今方程毅值解
上
机
实
验
报
(一)实验一
一.上机题目：
用线性元求解下列边值问题的数值解:
I, . 7T27T2. 7T小-I
-y H ----- y = —sin- x, 0 J 4 丿2 2 y(0)=0, y'(l) = o 二、实验程序： function S=bz x=fzero(@zfun r1); [t y]=ode45(@odefun r [0 1]z [0 x]); S.t=t; S.y=y; plot(t,y) xlabel(!x: z从0 ?直到1‘) ylabel (1 y1) title (，线性元求解边Fi问题的数值解?) function dy=odefun(x,y) dy=[0 0]?; dy (1) =y (2); dy (2) = (pi A2) /4*y (1) - ( (pi A2) /2) *sin (x*pi/2); function z=zfun(x); [t y]=ode45(@odefun z [0 1]z [0 x]); z=y(end)-0; 三、实验结果: 1.以步长h=0.05进行逐步运算，运行上面niatlab程序结果如下: ans t: [57x1 double： y: "57x2 doubled 2.在0 （二）实验二上机题目: 求解Hehnlioltz方程的边值问题： —Aw — k2u = 1,于Q = (0,1) *(0,1) u=o,于= {x = 050< y —=0,于I"2 = {0 S x S1, y = 0} U {x = 1,0 S y S 1} dn ' 其中k=l,5,10,15,20 五、实验程序： (采用冇限元方法，这里对[0,1]*[0,1]采用屮11的划分，n为偶数) n=10; a=zeros(n);f=zeros(n);b=zeros(l r n);U=zeros(n f1);u=zeros(n,1); for i=2:n a(i-l f i-l)=pi A2/(12*n)+n; a(i-l f i)= pi A2/(24*n)-n; a (i,i-1)= pi A2/(24*n)-n; for j=l:n if j==i-l a(i, j)=a(i,i-1)； else if j==i a (i-l f j-1)=2*a(i-1,i-1); else if j==i+l a(i, j)=a(i,i+1)； else a(i, j)=0; end end end end end a(n,n)=pi A2/(12*n)+n; for i=2:n f (i-1/ i) =4/pi*cos ( (i-1) *pi/2/n) -8*n/ (pi A2) *sin (i*pi/2/n) +8*n/ (pi A2) *s in ((i-1)*pi/2/n); end for i=l:n f (i z i)=-4/pi*cos(i*pi/2/n)+8*n/(pi A2)*sin(i*pi/2/n)-8*n/(pi A2)*sin((i -1)*pi/2/n); end %b(j)=f(i-1,j)+f(i, j) for i=l:(n-1) b(i)=f(i,i)+f (i r i+1); end b (n) =f (n f n); tic; n=20; can=20; s=zeros (22,10); h=l/n; st=l/(2*n A2); A=zeros( (n+1)A2# (n+1)A2); syms x y; for k=l:l:2*n A2 s(k,1)=k; q=fix(k/(2*n)); r=mod(k r (2*n)); if (r-=0) r=r; else r=2*n;q=q-l; end if (r<=n) s(k,2)=q*(n+1)+r; s(k,3)=q*(n+l)+r+l; s(k z4)=(q+l)*(n+l)+r+l; s(k z 5) = (r-1)*h; s(k,6)=q*h; s(k z7)=r*h; s(k z 8)=q*h; s(k,9)=r*h; s(k,10)=(q+l)*h; else s(k z 2)=q*(n+1)+r-n; s(k,3)=(q+1)*(n+1)+r-n+l; s(k z 4) = (q+1)*(n+1)+r-n; s(k,5)=(r-n-l)*h; s(k z 6)=q*h; s(k z7)=(r-n)*h; s(k,8)=(q+l)*h; s(k,9)=(r-n-l)*h; s(k,10)=(q+l)*h; end end d=zeros(3,3); B=zeros((n+1)A2f1); b=zeros (3,1); for k=l:l:2*n A2 L(l) = (l/(2*st))*( (s(k/7)*s(k r10)-s(k/9)*s(k r8)) + (s(k r8)-s(k f 10))*x+(s (k f9)-s(k r7))*y); L(2) = (l/(2*st) )*( (s(k/9)*s(k,6)-s(k z5)*s(k/10)) + (s(k r10)-s(k/6) )*x+(s (k,5)-s(k,9))*y); L(3) = (l/(2*st))*( (s(k/5)*s(k z8)-s(k z7)*s(k/6) ) + (s (k, 6)-s (k, 8) ) *x+(s (k ,7)-s(k,5))*y); for i=l:l:3 for j=i:3 d(i z j)=int(int(((((diff(L(i),x))*(diff(L(j),x)))+((diff(L(i),y))*(dif f (L(j) ,y) ) ) )-( (can A2) *L (i) *L (j ) ) ) , x, 0,1), y, 0,1); d(j, i)=d(i, j)； end end for i=l:l:3 for j=l:l:3 A(s(k, (i+l)),s(k, (j+1) ) ) =A(s (：c z (i+l)),s(k z (j+ 1) ) )+d (i, j ); end end for i=l:l:3 b (i) =int (int ( (L(i) ) ,x,0z 1) ,y,0,l); B(s(k,(i+1))r l)=B(s(k r (i+1)),l)+b(i)； end end M=zeros((n+1)A2r n A2); j=n A2; for i=(n A2+n):-l:l if ( (mod (i, (n+1)) ) -^=1) M(: r j ) =A (: r i)； else continue end end preanswer=M\B; answer=zeros((n+1)A2Z1); j=l; for i=l:l:(n A2+n) if ((mod(i f (n+1)))~=1) answer(i)=preanswer(j); j=j+l; else answer(i)=0; end end Z=zeros((n+1), (n+1)); for i=l:l:(n+1)A2 s=fix(i/(n+l))+l; r=mod(i,(n+1)); if(r==0) r=n+l; s=s-l; else end Z (r r s)=answer(i); end [X f Y]=meshgrid(l:-h:0,0:h:l); surf(X r Y r Z); toe; t=toc; K=a ; B=b? ; U=inv(K)*B for i=l:n u(i r l)=4/(pi A2)*sin(pi*i/n/2); end u e=U-u 六、实验结果： 程序中的变量can为问题中的k,为了方便比较，采用了画图的方式。 在n=10时，分别考察can=k=l,5,10,15,20时的情况。