五点差分格式资料.doc

五点差分格式《微分方程数值解》大作业(一)——椭圆型方程编程计算:采用五点差分格式求如下椭圆型方程其中、及边条件为:1.,且边条件如下:问题存在精确解为:2.,且边条件如下:问题存在精确解为:3.,且边条件如下:问题存在精确解为:.代码:主函数1,差分解functiong=fivepoints(x1,x2,y1,y2,M,N)%变步长法h=(x2-x1)/M;%横轴步长k=(y2-y21/N;%纵轴步长m=M-1;n=N-1;h1=h^2;r=h1/k^2;%五点中的上下两个点的系数t=2+2*r;%五点中的中心点的系数x=x1+(x2-x1)*(0:M)/M;%x,y向量表示横纵坐标y=y1+(y2-y1)*(0:N)/N;a=zeros(m*n,m*n);b=zeros(m*n,1);%初始化a,b矩阵,a为系数矩阵%内部的(m-2)*(n-2)个点fori=2:m-1forj=2:n-1a(i+(j-1)*m,:)=[zeros(1,i-1+(j-2)*m)-rzeros(1,m-2)-1t-1zeros(1,m-2)-rzeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1));endend%下边缘j=1;fori=2:m-1a(i+(j-1)*m,:)=[zeros(1,i-2)-1t-1zeros(1,m-2)-rzeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*bottom(x(i+1));end;%右边缘i=m;forj=2:n-1a(i+(j-1)*m,:)=[zeros(1,(j-1)*m-1)-rzeros(1,m-2)-1tzeros(1,m-1)-rzeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+right(y(j+1));end%上边缘j=n;fori=2:m-1a(i+(j-1)*m,:)=[zeros(1,i-1+(j-2)*m)-rzeros(1,m-2)-1t-1zeros(1,m-i-1)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*top(x(i+1));end%左边缘i=1;forj=2:n-1a(i+(j-1)*m,:)=[zeros(1,i-1+(j-2)*m)-rzeros(1,m-1)t-1zeros(1,m-2)-rzeros(1,(n-j)*m-i)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+left(y(j+1));end;%左下角的那个点i=1;j=1;a(1,:)=[t-1zeros(1,m-2)-rzeros(1,(n-1)*m-1)];b(1)=h1*f(x(2),y(2))+r*bottom(x(2))+left(y(2));%右下角的那个点i=m;j=1;a(i+(j-1)*m,:)=[zeros(1,m-2)-1tzeros(1,m-1)-rzeros(1,(n-2)*m)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*bottom(x(i+1))+right(y(j+1));%左上角的那个点i=1;j=n;a(i+(j-1)*m,:)=[zeros(1,(n-2)*m)-rzeros(1,m-1)t-1zeros(1,m-2)];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*top(x(i+1))+left(y(j+1));%右上角的那个点i=m;j=n;a(i+(j-1)*m,:)=[zeros(1,(n-1)*m-1)-rzeros(1,m-2)-1t];b(i+(j-1)*m)=h1*f(x(i+1),y(j+1))+r*top(x(i+1))+right(y(j+1));u=a\bab2,精确解:functiong=ni(x1,x2,y1,y2,M,N)m=M-1;n=N-1;x=x1+(x2-x1)*(0:M)/M;y=y1+(y2-y1)*(0:N)/N;fori=1:mforj=1:nu1(i+(j-1)*m)=f1(x(i+1),y(j+1))endend(1)辅助函数functiong=f(x,y)g=0;functiong=bottom(x)g=2*log(x);functiong=right(y)g=log(4+y^2);functiong=top(x)g=log(x^2+1);functiong=left(y)g=log(1+y^2);functiong=f1(x,y)g=log(x^2+y^2);运行fivepoints(1,2,0,1,4,4)u=