三、椭圆形方程算法和程序： 1、（a）用程序计算5*5的网格，确定9个未知数平p1、p2……p9的方程组，来求解矩形区域R={(x,y)|0≤x≤4,0≤y≤4}内的谐波函数u(x,y)的近似值。,边界为：
u(x,0)=10和u(x,4)=120,0 （b）用9*9的网格求解近似解。 2、用程序计算矩形区域R={(x,y)|0≤x≤1.5,0≤y≤1.5}内的谐波函数u(x,y)的近似值，h=0.15,边界为： u(x,0)=x4和u(x,4)=x4-13.5x2+5.0625,0 程序代码： function [u,w,x,y ] = dirich( f1,f2,f3,f4,a,b,h,tol,max1) %function：椭圆型方程的差分法 ?,f2,f3,f4：边界条件，字符型（string） %a，b：x，y上限值； %h：步长 n=fix(a/h)+1; m=fix(b/h)+1; ave=(a*(feval(f1,0)+feval(f2,0)) +b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b); U=ave*ones(n,m); %边界条件赋值 U(1,1:m)=feval(f3,0:h:(m-1)*h)'; U(n,1:m)=feval(f4,0:h:(m-1)*h)'; U(1:n,1)=feval(f1,0:h:(n-1)*h); U(1:n,m)=feval(f2,0:h:(n-1)*h); U(1,1)=(U(1,2)+U(2,1))/2; U(1,m)=(U(1,m-1)+U(2,m))/2; U(n,1)=(U(n-1,1)+U(n,2))/2; U(n,m)=(U(n-1,m)+U(n,m-1))/2; %SOR parameter（超松弛因子） w=4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2)); %计算 err=1; cnt=0; while((err>tol) & (cnt<=max1)) err=0; for j=2:m-1 for i=2:n-1 relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4*U(i,j))/4; U(i,j)=U(i,j)+relx; if (err<=abs(relx)) err=abs(relx); end end end cnt=cnt+1; end u=flipud(U'); end 问题1： （a） P1=54.2857 p2=41.4286 p3=36.4286 p4=75.7143 p5=65.0000 p6=54.2857 p7=93.5714 p8=88.5714 p9=75.7143 （b） 计算结果： 近似解图示 问题2：结果： 近似解图：