function L=Laplace_N(x0,y0,Lx,Ly,h)
I=Lx/h;J=Ly/h;
for(i=1:I+1)
    u(1,i)=feval('Bound_Con',x0+(i-1)*h,y0);
    u(J+1,i)=feval('Bound_Con',x0+(i-1)*h,y0+Ly);
end%计算边界值 先按行计算  边界条件:Boundary conditions
for(i=2:J)
    u(i,1)=feval('Bound_Con', x0,y0+(i-1)*h);
    u(i,I+1)=feval('Bound_Con',x0+Lx,y0+(i-1)*h);    %边界值计算 左右两列
end%求系数矩阵,除去边界点,内点一共有(I-1)^2个,根据系数矩阵的特点,我们采用单元数组的方式进行编程。
%先对模块矩阵A和B进行赋值
A=zeros(I-1);
for(i=1:I-1)
    for(j=1:I-1)
        if i==j
            A(i,j)=4;
        elseif abs(i-j)==1
            A(i,j)=-1;
        else
            A(i,j)=0;
        end
    end
end
B=-eye(I-1);
Z=zeros(I-1);%B赋值为单位矩阵
%下面对系数矩阵进行赋值利用单元数组的格式
C=cell(J-1,J-1);
for(i=1:J-1)
    for(j=1:J-1)
        if(i==j)
            C{i,j}=A;
        elseif(abs(i-j)==1)
            C{i,j}=B;
        else%if(i~=j&abs(i-j)~=1)
            C{i,j}=Z;
        end
    end
end
C1=cell2mat(C);%把单元数组转换为数值数组
for(i=1:(I-1)*(J-1))
    xi=(fix((i-1)/(I-1))+1)*h;
    yj=(i-(I-1)*fix((i-1)/(I-1)))*h;
        if(i==1)
            D(i)=u(1,2)+u(2,1)+h^2*feval('Fun',xi,yj);
        elseif(i==I-1)
            D(i)=u(1,I)+u(2,I+1)+h^2*feval('Fun',xi,yj);
        elseif(i==(I-1)*(J-2)+1)
            D(i)=u(J,1)+u(J+1,2)+h^2*feval('Fun',xi,yj);
        elseif(i==(I-1)*(J-1))
            D(i)=u(J+1,I)+u(J,I+1)+h^2*feval('Fun',xi,yj);%四个角点
        elseif(1<i&i<I-1)
            D(i)=u(1,i+1)+h^2*feval('Fun',xi,yj);
        elseif((J-2)*(I-1)+1<i&i<(I-1)*(J-1))
            D(i)=u(J+1,i-((I-1)*(J-1)-I))+h^2*feval('Fun',xi,yj)%上下两行的边界点值
        elseif(i-(I-1)*fix(i/(I-1))==0&i~=I-1&i~=(I-1)*(J-1))
            D(i)=u(fix(i/(I-1))+1,I+1)+h^2*feval('Fun',xi,yj);
        elseif(i-(I-1)*fix(i/(I-1))==1&i~=1&i~=(I-1)*(J-2))
            D(i)=u(fix(i/(I-1))+2,1)+h^2*feval('Fun',xi,yj);%左右两侧的边界点值
        else
            D(i)=0+h^2*feval('Fun',xi,yj);
        end
end%再对方程组右侧的常数项数组赋值
U=C1\D';%完成之后,进行矩阵除法运算
for(i=2:I)
    for(j=2:J)
        u(i,j)=U((i-2)*(I-1)+j-1);
end%此时,将U中的值赋回给u
surf(u);%绘图
end
end

点赞 ({{click_count}}) 收藏 (0)