兩點(diǎn)邊值微分方程的有限元分析報(bào)告

....重慶大學(xué)生實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)課程名稱偏微分方程數(shù)值解期末課程設(shè)計(jì) 開課實(shí)驗(yàn)室偏微分方程數(shù)值解學(xué)院數(shù)統(tǒng)學(xué)院年級2011專業(yè)班學(xué)生姓名學(xué)號 開課時(shí)間2013至2014學(xué)年第2學(xué)期總總成績教師簽名.數(shù)學(xué)與統(tǒng)計(jì)學(xué)院制驗(yàn)項(xiàng)目名稱指導(dǎo)教師偏微分方程期末課程設(shè)計(jì)成績實(shí)驗(yàn)項(xiàng)目類型驗(yàn)證演示綜合設(shè)計(jì)其他。。實(shí)驗(yàn)容。 ....三、實(shí)驗(yàn)原理、方法(算法)、步驟SobolevH(a，b)的子空間V:V={|∈H1(a，b)，v(a)=0}令....…有限元方如下性質(zhì)顯然，基函數(shù)線性無關(guān)，對任意的一個(gè)e可表示為因此，求e，使得相當(dāng)于：求相當(dāng)于：求.......xnnxh2h2n1x3x33xx2+3x2x3x+3xnnnnnnnn1i1ixh2h2iiiii1i1x3iiiii1ih6h2i=2Kna(Q,Q)=jxi+11(xi+1x)(xxi)dxi=1Kn-1ii+1xh2h2i=i+1iii+1ii+11x3x3=i+1iii+1ii+1h6h2.ixhxh=i1iiiii=i1iiiii1i1i+1i+1niii=1實(shí)驗(yàn)二;入2u+2(1+入2)u入2uuu=0i+1,ji,ji1,ji,j+1i,j1(1)x=0:0,=0,forj=1,...,M,(2)x=1:+1,=0,forj=1,...,M,(3)y=0:,0=0,fori=1,...,N,(4)y=1:,+1=(),fori=1,...,N.可以改寫為矩陣形式為：....行了，由此可以看出，系數(shù)矩陣A隨著網(wǎng)格數(shù)的增加急劇的變大以至于電腦存存不的矩陣，而令bx+ax+bx+cx=y(i)i1ii+1N+ibx+ax+cx=y(i)i1iN+icx+ax+bx+cx=y(i)iNii+1i+N當(dāng)....當(dāng)當(dāng)cx+bx+ax+bx=y(i)iNi1ii+1cx+bx+ax=y(i)iNi1i因此這樣就解決了矩陣A存儲不下的尷尬情況，大大的減少了程序所占存的大四．實(shí)驗(yàn)環(huán)境(所用軟件、硬件等)及實(shí)驗(yàn)數(shù)據(jù)文件MATLAB14afunctionmyfun1(n)h=1/n;x=0:h:1;A=zeros(n,n);b=zeros(n,1);fori=2:n%生成剛性矩陣A(i-1,i-1)=-(x(i-1)^3+3*x(i-1)^2*x(i)-3*x(i-1)*x(i)^2-2*x(i)^3+3*x(i)^2*x(i+1)-3*x(i)*x(i+1)^2+x(i+1)^3-3*x(i-1)+3*x(i+1))/(3*h^2);A(i-1,i-1)=A(1,1);A(n,n)=-1/h+1/h^2*(x(n+1)^3/3-x(n)^3/3-x(n)*x(n+1)^2+x(n)^2*x(n+1));A(n,n)=A(1,1)/2;fori=2:nA(i-1,i)=1/h-1/h^2*(x(i+1)^3/6-x(i)^3/6+x(i)^2*x(i+1)/2-x(i)*x(i+1)^2/2);A(i,i-1)=A(i-1,i);fori=2:n%生成右端常數(shù)項(xiàng)b(i-1)=1/(2*h*pi^2)*(-x(i-1)*cos(2*pi*x(i))*pi+2*cos(2*pi*x(i))*pi*x(i)-x(i+1)*cos(2*pi*x(i))*pi+sin(pi*x(i-1))*cos(pi*x(i-1))+sin(pi*x(i+1))*cos(pi*x(i+1))-sin(2*pi*x(i)));b(n)=1/(4*h*pi^2)*(-2*x(n)*cos(2*pi*x(n+1))*pi+2*cos(2*pi*x(n+1))*pi*x(n+1)+2*sin(pi*x(n))*cos(pi*x(n))-sin(2*pi*x(n+1)))-2.03*pi/(1+4*pi^2);y=myconjgrad(A,b);%調(diào)用求解方程組的算法u1=zeros(n+1,1);....u1(2:n+1)=y;x=0:0.01:1;y=sin(2*pi*x)/(1+4*pi^2);legend('有限元','精確解');.....functionmyElliptic1(N)h=1/(N+1);k=1/(N+1);l=k/h;x=0:h:1;y=0:k:1;g=sin(5*pi*x);a=zeros(N,N);%生成系數(shù)矩陣b=ones(N,1)*(2*(1+l^2));c=ones(N-1,1)*(-l^2);a=diag(b)+diag(c,1)+diag(c,-1);A=zeros(N^2,N^2);fori=1:NA((i-1)*N+1:i*N,(i-1)*N+1:i*N)=a;A=A+diag(-ones(N*(N-1),1),N)+diag(-ones(N*(N-1),1),-N);B=zeros(N*N,1);%生成右端常向量B(N*(N-1)+1:N*N)=g(2:N+1);v=myconjgrad(A,B);u=zeros(N,N);u(:)=v;m=zeros(N+2,N+2);m(:,N+2)=g;m(2:N+1,2:N+1)=u;surf(x(N+2:-1:1),y,m);holdon;xlabel('x-axis');xlabel('y-axis');zlabel('Siolution');Elliptic1(120)A=A+diag(-ones(N*(N-1),1),N)+diag(-ones(N*(N-1),1),-N);functionmyElliptic2(N)h=1/(N+1);k=1/(N+1);l=k/h;x=0:h:1;y=0:k:1;....g=sin(5*pi*x);%只存儲系數(shù)矩陣中的三個(gè)元素a1=(2*(1+l^2));b1=(-l^2);c%只存儲系數(shù)矩陣中的三個(gè)元素B=zeros(N*N,1);%調(diào)用共軛梯度法計(jì)算線性方程組的解B(N*(N-1)+1:N*N)=g(2:N+1)%調(diào)用共軛梯度法計(jì)算線性方程組的解v=myconjgrad1(N,a1,b1,c1,B);u=zeros(N,N);u(:)=v;m=zeros(N+2,N+2);m(:,N+2)=g;%作圖m(2:N+1,2:N+1)=u;%作圖surf(x(N+2:-1:1),y,m);holdon;xlabel('x-axis');xlabel('y-axis');zlabel('Siolution');functionx=myconjgrad1(N,a1,b1,c1,B)a2=a1;b2=b1;c2=c1;N2=N;b=B;%設(shè)定誤差限%調(diào)用壓縮存貯后的%設(shè)定誤差限%調(diào)用壓縮存貯后的A*b的算法%設(shè)定循環(huán)次數(shù)x=b;r=b-fun(N2,a2,b2,c2,x);p=r;fork=1:(numel(b))^2;z=myfun(N2,a2,b2,c2,p);alpha=(r'*r)/(p'*z);x=x+alpha*p;s=r'*r;r=r-alpha*z;if(norm(r)<tol)return;B=(r'*r)/s;p=r+B*p;functionv=myfun(N,a,b,c,x)v=zeros(N*N,1);v(1)=a*x(1)+b*x(2)+c*x(N+1);fori=2:N-1;v(i)=b*x(i-1)+a*x(i)+b*x(i+1)+c*x(i+N);v(N)=b*x(N-1)+a*x(N)+c*x(N+N);v(N+1)=c*x(1)+a*x(N+1)+b*x(N+2)+c*x(N+1+N);fori=N+2:N*N-N-1;v(i)=c*x(i-N)+b*x(i-1)+a*x(i)+b*