[已答复] poisson方程matlab程序!求解区域是长方形额。

[复制链接]
四季之梦 发表于 2014-6-13 21:45:02
就是这个方程。。然后之前有写过区间是【0,1】*【0,1】的。。但是这个区间不是正方形的构造出的矩阵就不是(n-1)^2*(n-1)^2而是((m-1)*(n-1))*((m-1)*(n-1))了。。这样就不知道左端的 矩阵长什么样子了。。有大神知道吗??急求啊。。function [e] = problem5(M)h=1/M;%(这个之前我写的另一个是正方形区域的程序。。求帮忙长方形区域的程序)
for i=1:M-1%剖分,得内点
    x(i)=i*h;
    y(i)=i*h;
end
A=zeros(M-1,M-1);%初始
B=zeros((M-1)^2,(M-1)^2);
E=-1*eye(M-1);%生成小块对角矩阵
for i=1:M-1%生成小块三对角矩阵
    for j=1:M-1
if i==j
    A(i,j)=4;
else if i+1==j|i==j+1
       A(i,j)=-1;
    end
end
    end
end
for i=1:M-1
for j=1:M-1
        f(i,j)=sin(pi*x(i))*sin(pi*y(j))+sin(pi*x(i))*sin(2*pi*y(j)); %生成f矩阵,
        u(i,j)=(2*pi^2)^(-1)*sin(pi*x(i))*sin(pi*y(j))+(5*pi^2)^(-1)*sin(pi*x(i))*sin(2*pi*y(j));%生成真实解矩阵
    end
end
for k=1:M-1%将小块矩阵A放到B里面
    for i=1+(k-1)*(M-1):k*(M-1)
        b(i,1)=f(i-(k-1)*(M-1),k);%把f弄成长向量
        for j=1+(k-1)*(M-1):k*(M-1)
            B(i,j)=A(i-(k-1)*(M-1),j-(k-1)*(M-1));
        end
    end
end
for k=1:M-2%将小块矩阵E放到B里面
    for i=1+(k-1)*(M-1):k*(M-1)
        for j=i+(M-1):(k+1)*(M-1)
            B(i,j)=E(i-(k-1)*(M-1),j-k*(M-1));
        end
    end
end
for i=1+(M-2)*(M-1):(M-1)^2
     for j=1+(M-3)*(M-1):(M-2)*(M-1)
            B(i,j)=E(i-(M-2)*(M-1),j-(M-3)*(M-1));
     end
end
U=h^2*b'/B;%解方程组
for k=1:M-1%将U弄成矩阵
    for i=1+(k-1)*(M-1):k*(M-1)
        UU(i-(k-1)*(M-1),k)=U(1,i);
    end
end
e=u(5,5)-UU(5,5);%在点(0.5,0.5)处的误差
u00=zeros(M+1,1);u10=zeros(1,M-1);%边界,把边界值加上去
uu=[u00 [u10;u;u10] u00];
UU=[u00 [u10;UU;u10] u00];
xx=[0 x 1];
yy=[0 y 1];
mesh(xx,yy,uu)%画图
figure
mesh(xx,yy,UU)%画图

end



QQ图片20140613214236.jpg

1 条回复


四季之梦 发表于 2014-6-13 21:50:09
[e] = problem5(M);
h=1/M;前两句打上去的时候不小心弄错了。。。我给的程序是【0,1】*【0,1】区间的。。我现在要做的事【1,2】*【0,2pi】的。。额。。方程在图片里。。用的五点差分格式。。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表