偏微分方程數(shù)值解實(shí)驗(yàn)報(bào)告

--上上機(jī)實(shí)驗(yàn)報(bào)告----(一)實(shí)驗(yàn)一-?′′+?2?=?2sin??，0<x<1422y(0)=0,?′(1)=0functionS=bzx=fzero(＠zfun,1)；[ty]＝ode45(@odefun,[01],[0x]);S.t=t;S.y＝y(tǒng);plot(t，y)xlabel('x:′從0一直到1')ylabel('y')title('線性元求解邊值問題的數(shù)值解')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;--(二)實(shí)驗(yàn)二求解Helmholtz方程的邊值問題：1?u=0,于={0x1,y=0}{x=1,0y1}?n?n----(采用有限元方法，這里對[0，1]*[0,1]采用n*n的劃分,n為偶數(shù))n=10;a=zeros(n);f=zeros(n);b=zeros(1,n);U=zeros(n,1);u=zeros(n，1)；fori＝2:na(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；orjnifj==i-1a(i，j)=a(i,i-1);elseifj＝=ia(i-1,j－1)=2*a(i-1,i－1);elseifj==i+1a(i，j)=a(i,i+1);elsea(i,j)＝0;endendendenda(n，n)=pi^2/(12*n)+n;fori＝2:nf(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);endfori=1:nf(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);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^2s(k，1)＝k;q=fix(k/(2＊n));r=mod(k,(2*n));if(r~=0)r=r;elser=2*n；q=q-1;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；elses(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；endendd=zeros(3,3);B＝zeros((n+1)^2,1)；b=zeros(3,1);fork=1:1:2*n^2L(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:3forj=i:3d(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)；endfori=1:1:3forj=1:1：3A(s(k,(i+1)),s(k,(j+1)))=A(s(k，(i+1)),s(k,(j+1)))＋d(i，j)；endendfori=1:1:3b(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);endendM=zeros((n+1)^2,n^2);j＝n^2；fori=(n^2+n)：－1：1if((mod(i,(n+1)))~＝1)M(：,j)＝A(:,i);j=j-1；elsecontinueendpreanswer=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;endendZ=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;elseZ(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)*Bfori=1:nu(