我在学习MATLAB的使用,想用有限差分的方法求解二阶偏微分方程,请高手指点!

时间:2022-02-26 19:31:22
我已推导出有限差分公式,可不知如何用MATLAB求解,请指点,多谢!最好能给出源代码。方程式为什么粘贴不上?

12 个解决方案

#1


没看到具体方程,不过这个应该对你有用。

用五点差分格式解拉普拉斯方程


function u = peEllip5(nx,minx,maxx,ny,miny,maxy)
format long;
hx = (maxx-minx)/(nx-1);
hy = (maxy-miny)/(ny-1);
u0 = zeros(nx,ny);
for j=1:ny
    u0(j,1) = EllIni2Uxl(minx,miny+(j-1)*hy);
    u0(j,nx) = EllIni2Uxr(maxx,miny+(j-1)*hy);
end
for j=1:nx
    u0(1,j) = EllIni2Uyl(minx+(j-1)*hx,miny);
    u0(ny,j) = EllIni2Uyr(minx+(j-1)*hx,maxy);
end 

A = -4*eye((nx-2)*(ny-2),(nx-2)*(ny-2));
b = zeros((nx-2)*(ny-2),1);
for i=1:(nx-2)*(ny-2)
    if mod(i,nx-2) == 1
        if i==1
            A(1,2) = 1;
            A(1,nx-1) = 1;
            b(1) = - u0(1,2) - u0(2,1);
        else
            if i == (ny-3)*(nx-2)+1
                A(i,i+1) = 1;
                A(i,i-nx+2) = 1;
                b(i) = - u0(ny-1,1) - u0(ny,2);
            else
                A(i,i+1) = 1;
                A(i,i-nx+2) = 1;
                A(i,i+nx-2) = 1;
                b(i) = - u0(floor(i/(nx-2))+2,1);
            end
        end
    else
        if mod(i,nx-2) == 0
            if i == nx-2
                A(i,i-1) = 1;
                A(i,i+nx-2) = 1;
                b(i) = - u0(1,nx-1) - u0(2,nx);
            else
                if i == (ny-2)*(nx-2)
                    A(i,i-1) = 1;
                    A(i,i-nx+2) = 1;
                    b(i) = - u0(ny-1,nx) - u0(ny,nx-1);
                else
                    A(i,i-1) = 1;
                    A(i,i-nx+2) = 1;
                    A(i,i+nx-2) = 1;
                    b(i) = - u0(floor(i/(nx-2))+1,nx);
                end
            end
        else
            if i>1 && i< nx-2
                A(i,i-1) = 1;
                A(i,i+nx-2) = 1;
                A(i,i+1) = 1;
                b(i) = - u0(1,i+1);
            else
                if i > (ny-3)*(nx-2) && i < (ny-2)*(nx-2)
                    A(i,i-1) = 1;
                    A(i,i-nx+2) = 1;
                    A(i,i+1) = 1;
                    b(i) = - u0(ny,mod(i,(nx-2))+1);
                else
                    A(i,i-1) = 1;
                    A(i,i+1) = 1;
                    A(i,i+nx-2) = 1;
                    A(i,i-nx+2) = 1;
                end
            end
        end
    end
end
ul = A\b;
for i=1:(ny-2)
    for j=1:(nx-2)
        u(i,j) = ul((i-1)*(nx-2)+j);
    end
end
format short;    

#2


用工字型差分格式解拉普拉斯方程


function u = peEllip5m(nx,minx,maxx,ny,miny,maxy)
format long;
hx = (maxx-minx)/(nx-1);
hy = (maxy-miny)/(ny-1);
u0 = zeros(nx,ny);
for j=1:ny
    u0(j,1) = EllIni2Uxl(minx,miny+(j-1)*hy);
    u0(j,nx) = EllIni2Uxr(maxx,miny+(j-1)*hy);
end
for j=1:nx
    u0(1,j) = EllIni2Uyl(minx+(j-1)*hx,miny);
    u0(ny,j) = EllIni2Uyr(minx+(j-1)*hx,maxy);
end 

A = -4*eye((nx-2)*(ny-2),(nx-2)*(ny-2));
b = zeros((nx-2)*(ny-2),1);
for i=1:(nx-2)*(ny-2)
    if mod(i,nx-2) == 1
        if i==1
            A(1,nx) = 1; 
            b(1) = - u0(1,1) - u0(3,1) - u0(1,3);
        else
            if i == (ny-3)*(nx-2)+1
                A(i,i-nx+1) = 1;
                b(i) = - u0(ny,1) - u0(ny,3) - u0(ny-2,1);
            else
                A(i,i-nx+3) = 1;
                A(i,i+nx-1) = 1;
                b(i) = - u0(floor(i/(nx-2))+1,1) - u0(floor(i/(nx-2))+3,1);
            end
        end
    else
        if mod(i,nx-2) == 0
            if i == nx-2
                A(i,i+nx-1) = 1;
                b(i) = - u0(1,nx-2) - u0(1,nx) - u0(3,nx);
            else
                if i == (ny-2)*(nx-2)
                    A(i,i-nx+1) = 1;
                    b(i) = - u0(ny,nx) - u0(ny,nx-2) - u0(ny-2,nx);
                else
                    A(i,i-nx+1) = 1;
                    A(i,i+nx-3) = 1;
                    b(i) = - u0(floor(i/(nx-2)),nx) - u0(floor(i/(nx-2))+2,nx);
                end
            end
        else
            if i>1 && i< nx-2
                A(i,i+nx-1) = 1;
                A(i,i+nx-3) = 1;
                b(i) = - u0(1,i) - u0(1,i+2);
            else
                if i > (ny-3)*(nx-2) && i < (ny-2)*(nx-2)
                    A(i,i-nx+3) = 1;
                    A(i,i-nx+1) = 1;
                    b(i) = - u0(ny,mod(i,(nx-2))) -u0(ny,mod(i,(nx-2))+2);
                else
                    A(i,i-nx+3) = 1;
                    A(i,i+nx-1) = 1;
                    A(i,i+nx-3) = 1;
                    A(i,i-nx+1) = 1;
                end
            end
        end
    end
end
ul = A\b;
for i=1:(ny-2)
    for j=1:(nx-2)
        u(i,j) = ul((i-1)*(nx-2)+j);
    end
end
format short;    

#3


mark

#4


太好了!非常感谢!我可以把你加为好友吗?

#5


强,顶一个

#6


引用 4 楼 JENNY1999 的回复:
太好了!非常感谢!我可以把你加为好友吗?


当然可以~~

#7


我已加你为好友,并有留言,请关注,谢谢!

#8


你发给我的程序运行时总是出现问题:

??? function u = peEllip5(nx,minx,maxx,ny,miny,maxy)
    |
Error: Function definitions are not permitted at the prompt or in scripts.

你能否帮我解决?谢谢!

#9


matlab有函数的

#10


up

#11


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#12


能教我这个新手怎么加大家好友么?我的qq964707133,想加这里的好友,帮忙呀,谢谢大家了

#1


没看到具体方程,不过这个应该对你有用。

用五点差分格式解拉普拉斯方程


function u = peEllip5(nx,minx,maxx,ny,miny,maxy)
format long;
hx = (maxx-minx)/(nx-1);
hy = (maxy-miny)/(ny-1);
u0 = zeros(nx,ny);
for j=1:ny
    u0(j,1) = EllIni2Uxl(minx,miny+(j-1)*hy);
    u0(j,nx) = EllIni2Uxr(maxx,miny+(j-1)*hy);
end
for j=1:nx
    u0(1,j) = EllIni2Uyl(minx+(j-1)*hx,miny);
    u0(ny,j) = EllIni2Uyr(minx+(j-1)*hx,maxy);
end 

A = -4*eye((nx-2)*(ny-2),(nx-2)*(ny-2));
b = zeros((nx-2)*(ny-2),1);
for i=1:(nx-2)*(ny-2)
    if mod(i,nx-2) == 1
        if i==1
            A(1,2) = 1;
            A(1,nx-1) = 1;
            b(1) = - u0(1,2) - u0(2,1);
        else
            if i == (ny-3)*(nx-2)+1
                A(i,i+1) = 1;
                A(i,i-nx+2) = 1;
                b(i) = - u0(ny-1,1) - u0(ny,2);
            else
                A(i,i+1) = 1;
                A(i,i-nx+2) = 1;
                A(i,i+nx-2) = 1;
                b(i) = - u0(floor(i/(nx-2))+2,1);
            end
        end
    else
        if mod(i,nx-2) == 0
            if i == nx-2
                A(i,i-1) = 1;
                A(i,i+nx-2) = 1;
                b(i) = - u0(1,nx-1) - u0(2,nx);
            else
                if i == (ny-2)*(nx-2)
                    A(i,i-1) = 1;
                    A(i,i-nx+2) = 1;
                    b(i) = - u0(ny-1,nx) - u0(ny,nx-1);
                else
                    A(i,i-1) = 1;
                    A(i,i-nx+2) = 1;
                    A(i,i+nx-2) = 1;
                    b(i) = - u0(floor(i/(nx-2))+1,nx);
                end
            end
        else
            if i>1 && i< nx-2
                A(i,i-1) = 1;
                A(i,i+nx-2) = 1;
                A(i,i+1) = 1;
                b(i) = - u0(1,i+1);
            else
                if i > (ny-3)*(nx-2) && i < (ny-2)*(nx-2)
                    A(i,i-1) = 1;
                    A(i,i-nx+2) = 1;
                    A(i,i+1) = 1;
                    b(i) = - u0(ny,mod(i,(nx-2))+1);
                else
                    A(i,i-1) = 1;
                    A(i,i+1) = 1;
                    A(i,i+nx-2) = 1;
                    A(i,i-nx+2) = 1;
                end
            end
        end
    end
end
ul = A\b;
for i=1:(ny-2)
    for j=1:(nx-2)
        u(i,j) = ul((i-1)*(nx-2)+j);
    end
end
format short;    

#2


用工字型差分格式解拉普拉斯方程


function u = peEllip5m(nx,minx,maxx,ny,miny,maxy)
format long;
hx = (maxx-minx)/(nx-1);
hy = (maxy-miny)/(ny-1);
u0 = zeros(nx,ny);
for j=1:ny
    u0(j,1) = EllIni2Uxl(minx,miny+(j-1)*hy);
    u0(j,nx) = EllIni2Uxr(maxx,miny+(j-1)*hy);
end
for j=1:nx
    u0(1,j) = EllIni2Uyl(minx+(j-1)*hx,miny);
    u0(ny,j) = EllIni2Uyr(minx+(j-1)*hx,maxy);
end 

A = -4*eye((nx-2)*(ny-2),(nx-2)*(ny-2));
b = zeros((nx-2)*(ny-2),1);
for i=1:(nx-2)*(ny-2)
    if mod(i,nx-2) == 1
        if i==1
            A(1,nx) = 1; 
            b(1) = - u0(1,1) - u0(3,1) - u0(1,3);
        else
            if i == (ny-3)*(nx-2)+1
                A(i,i-nx+1) = 1;
                b(i) = - u0(ny,1) - u0(ny,3) - u0(ny-2,1);
            else
                A(i,i-nx+3) = 1;
                A(i,i+nx-1) = 1;
                b(i) = - u0(floor(i/(nx-2))+1,1) - u0(floor(i/(nx-2))+3,1);
            end
        end
    else
        if mod(i,nx-2) == 0
            if i == nx-2
                A(i,i+nx-1) = 1;
                b(i) = - u0(1,nx-2) - u0(1,nx) - u0(3,nx);
            else
                if i == (ny-2)*(nx-2)
                    A(i,i-nx+1) = 1;
                    b(i) = - u0(ny,nx) - u0(ny,nx-2) - u0(ny-2,nx);
                else
                    A(i,i-nx+1) = 1;
                    A(i,i+nx-3) = 1;
                    b(i) = - u0(floor(i/(nx-2)),nx) - u0(floor(i/(nx-2))+2,nx);
                end
            end
        else
            if i>1 && i< nx-2
                A(i,i+nx-1) = 1;
                A(i,i+nx-3) = 1;
                b(i) = - u0(1,i) - u0(1,i+2);
            else
                if i > (ny-3)*(nx-2) && i < (ny-2)*(nx-2)
                    A(i,i-nx+3) = 1;
                    A(i,i-nx+1) = 1;
                    b(i) = - u0(ny,mod(i,(nx-2))) -u0(ny,mod(i,(nx-2))+2);
                else
                    A(i,i-nx+3) = 1;
                    A(i,i+nx-1) = 1;
                    A(i,i+nx-3) = 1;
                    A(i,i-nx+1) = 1;
                end
            end
        end
    end
end
ul = A\b;
for i=1:(ny-2)
    for j=1:(nx-2)
        u(i,j) = ul((i-1)*(nx-2)+j);
    end
end
format short;    

#3


mark

#4


太好了!非常感谢!我可以把你加为好友吗?

#5


强,顶一个

#6


引用 4 楼 JENNY1999 的回复:
太好了!非常感谢!我可以把你加为好友吗?


当然可以~~

#7


我已加你为好友,并有留言,请关注,谢谢!

#8


你发给我的程序运行时总是出现问题:

??? function u = peEllip5(nx,minx,maxx,ny,miny,maxy)
    |
Error: Function definitions are not permitted at the prompt or in scripts.

你能否帮我解决?谢谢!

#9


matlab有函数的

#10


up

#11


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#12


能教我这个新手怎么加大家好友么?我的qq964707133,想加这里的好友,帮忙呀,谢谢大家了