微分方程数值解实验MATLAB——线性椭圆型方程有限差分法

 拉普拉斯方程第一边值问题，五点差分法

核心思想就是剖分，将剖分得到的节点编号成一个向量，然后利用差分方程列方程组，最后求解即可。可能有不对的地方，请谅解，也期待有心人的改正。

hx = 1/3;%步长
hy = 1/3;
nx = 4;%每个方向的点数
ny = 4;
A = zeros(nx*ny, nx*ny);%系数矩阵A和B，所有点看做未知x
B = A;
U = zeros(nx*ny,1);%初值
for i = 0:ny-1
    for j = 0:nx-1
        I = nx*i + j + 1;%对点ij编号，转化成向量
        if mod(i,ny-1) == 0 || mod(j,nx-1) == 0%是边值点
            %边值条件得到方程组 Bx = U
            B(I,I) = 1;
            U(I,1) = bian(j*hx,i*hy);
        else 
            %差分方程得到的方程组 Ax = 0
            A(I,I) = -1;%点ij

            i = i - 1;%点i-1,j
            J = nx*i + j + 1;
            A(I,J) = 1/4;

            i = i + 2;%点i+1,j
            J = nx*i + j + 1;
            A(I,J) = 1/4;

            i = i - 1;
            j = j - 1;%点i,j-1
            J = nx*i + j + 1;
            A(I,J) = 1/4;

            j = j + 2;%点i,j+1
            J = nx*i + j + 1;
            A(I,J) = 1/4;            

        end
    end
end

x = (B-A) \ U;%求解方程组，得到所有点的值x
C = [];%最终解
for i = 0:ny-1%把向量x转换成矩阵
    C(i+1,:) = x(i*nx+1:i*nx+4)';
end

function u = bian(x,y)%边值条件
u = log((1+x)^2 + y^2);
end

运行结果

左上角为坐标原点，横轴为x，纵轴为y

第二个问题

hx = pi/5;
hy = 1/5;
nx = 6;
ny = 6;
A = zeros(nx*ny, nx*ny);
B = A;
U = zeros(nx*ny,1);
C = U;
tu = [];%精确解
for i = 0:ny-1
    for j = 0:nx-1
        tu(i+1,j+1) = trueU(j*hx,i*hy);
        I = nx*i + j + 1;
        if mod(i,ny-1) == 0 %y边界
            B(I,I) = 1; %是边值点
            U(I,1) = bian(j*hx,i*hy);
        elseif j == 0 %x边界
            B(I,I) = -1;
            B(I,I+1) = 1;
        elseif j == nx-1 %x边界
            B(I,I) = -1;
            B(I,I-1) = 1;
        else
            A(I,I) = 2*(hx^2 + hy^2);
            C(I,1) = f(j*hx,i*hy)*hx^2*hy^2;

            i = i - 1;
            J = nx*i + j + 1;
            A(I,J) = -hx^2;

            i = i + 2;
            J = nx*i + j + 1;
            A(I,J) = -hx^2;

            i = i - 1;
            j = j - 1;
            J = nx*i + j + 1;
            A(I,J) = -hy^2;

            j = j + 2;
            J = nx*i + j + 1;
            A(I,J) = -hy^2;            

        end
    end
end

x = (B-A)\(U - C);
u = [];%数值解
for i = 0:ny-1
    u(i+1,:) = x(i*nx+1:i*nx+nx)';
end

x = 0:hx:pi;
y = 0:hy:1;
[x,y] = meshgrid(x,y);
mesh(x,y,tu-u);

function a = f(x,y)
a = cos(3*x)*sin(pi*y);
end

function u = bian(x,y)
u = 0;
end

function u = trueU(x,y)
u = cos(3*x)*sin(pi*y)/(9+pi^2);
end

运行结果

九点差分格式

#九点，正方形
hx = 2/10;
hy = 2/10;
nx = 11;
ny = 11;
A = zeros(nx*ny, nx*ny);
B = A;
U = zeros(nx*ny,1);
C = U;
a = (hx^2 + hy^2)/(hx^2 * hy^2);
for i = 0:ny-1
    for j = 0:nx-1
        I = nx*i + j + 1;
        if mod(i,ny-1) == 0 || mod(j,nx-1) == 0
            B(I,I) = 1;
            U(I,1) = bian(j*hx-1,i*hy-1);
        else
            A(I,I) = 2/hx^2 + 2/hy^2 - a/3;
            C(I,1) = f(j*hx-1,i*hy-1);

            i = i - 1;
            j = j - 1;
            J = nx*i + j + 1;
            A(I,J) = -a/12;

            j = j + 1;
            J = nx*i + j + 1;
            A(I,J) = -1/hx^2 + a/6;

            j = j + 1;
            J = nx*i + j + 1;
            A(I,J) = -a/12;

            i = i + 1;
            j = j - 2;
            J = nx*i + j + 1;
            A(I,J) = -1/hy^2 + a/6;

            j = j + 2;
            J = nx*i + j + 1;
            A(I,J) = -1/hy^2 + a/6;            

            i = i + 1;
            j = j - 2;
            J = nx*i + j + 1;
            A(I,J) = -a/12; 

            j = j + 1;
            J = nx*i + j + 1;
            A(I,J) = -1/hx^2 + a/6;

            j = j + 1;
            J = nx*i + j + 1;
            A(I,J) = -a/12;

            i = i - 1;
        end
    end
end

x = (B-A)\(U - C);
u = [];
for i = 0:ny-1
    u(i+1,:) = x(i*nx+1:i*nx+nx)';
end

function a = f(x,y)
a = 16;
end

function u = bian(x,y)
u = 0;
end

运行结果