1.稳态温度场的分布(拉普拉斯方程第一边值问题数值解)
已有 665 次阅读2010-10-13 01:21|个人分类:课程实验|系统分类:科研笔记|关键词:laplace equation, numerical resolve
需要上机练习编程：差分法解拉普拉斯方程的第一边值问题。
自己编制的程序如下：
文件名：Lap-Eq Numerical answer.m
clc;clear;
tic
N=50
%划分的网格数======================
for m=1:N n=1:N-1;
u(m,n)=0;
u(m,N)=sin((m-1)*pi/(N-1));
end
%定义边界条件=======================
delta=ones(N,N);
while delta>1e-6
for m=2:N-1 n=2:N-1;
a(m,n)=u(m,n);
u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4;
delta(m,n)=abs(u(m,n)-a(m,n))/u(m,n);
end
end
X=1:N;Y=1:N;
mesh(X,Y,u(X,Y))
toc
所用的计算时间为Elapsed time is 3.672000 seconds.
1.考虑程序中的循环控制条件“while delta>=10e-6”的意义。
经过单步调试，得知这个表达式只是对最后一个delta进行比较，而不是所有的delta，因此并不满足计算条件。
结果是错误的。
要求每个计算点的delta都要<10e-6，因此需要该在程序。
clc;clear;
tic
N=50
%划分的网格数======================
for m=1:N n=1:N-1;
u(m,n)=0;
u(m,N)=sin((m-1)*pi/(N-1));
end
%定义边界条件=======================
delta=ones(N,N);
for m=2:N-1 n=2:N-1;
while delta(m,n)>1e-4
for m=2:N-1 n=2:N-1;
a(m,n)=u(m,n);
u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4;
delta(m,n)=abs(u(m,n)-a(m,n))/u(m,n);
end
end
end
X=1:N;Y=1:N;
mesh(X,Y,u(X,Y))
delta
delta>1e-6
toc
这样一来，所有的点都满足了。
但是这种算法做了太多的冗余计算。
对每个点的delta分别调到误差范围，所作的计算次数太多太多了。
从时间可看出
Elapsed time is 144.610000 seconds.需要进行算法改进，考虑每次计算结束得到一个delta的矩阵，只要矩阵中的最大者满足误差范围则所有的点都满足了，因此改为：
clc;clear;
tic
N=50
%划分的网格数======================
for m=1:N n=1:N-1;
u(m,n)=0;
u(m,N)=sin((m-1)*pi/(N-1));
end
%定义边界条件=======================
delta=zeros(N,N);
maxd=1;
whilemaxd>1e-4
for m=2:N-1 n=2:N-1;
a(m,n)=u(m,n);
u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4;
delta(m,n)=abs(u(m,n)-a(m,n))/u(m,n);
end
maxd=max(delta(:));
end
X=1:N;Y=1:N;
mesh(X,Y,u(X,Y))
maxd<1e-4
toc
现在就完全解决了上述问题了。
Elapsed time is 1.954000 seconds.
若误差要求改为1e-5，则运行时间为Elapsed time is 3.797000 seconds.
若误差要求改为1e-6，则运行时间为Elapsed time is 5.687000 seconds.
若划分的网格节点数N=500，tolerance=1e-5=> Elapsed time is 3794.234000 seconds.
几点需要说明的：
1) 对于二维或多维矩阵，找其最大值的表达式为max(A(:))，A代表矩阵名称。
A(:)代表矩阵A的所有元素以单序号方式引用。
这样找到的最大值才是一个数值。
若单纯的使用max(A)则对于二维矩阵会得到一个行矩阵，对应于A中每列的最
大值。
2) maxd<1e-4的作用，是检验是否所有的delta都已经满足误差要求了。
若满足，该式子的返回值为1，即为真。
另外，改变delta的存储情况也可减少存储空间，加快计算。
以下是不用矩阵存储delta，因为我们不需要知道每个delta值的表现形式，因此可以对每个delta 进行比较只用一个值来存储它。
程序如下：
%椭圆型方程的数值计算典型例题，Laplace方程的第一边值问题。
clc;clear;clf;
tic
N=5 %划分的网格节点数
tol=1e-5 %差分误差tolerance要求
%计算精度控制参量======================
for m=1:N n=1:N-1;
u(m,n)=0;
u(m,N)=sin((m-1)*pi/(N-1));
end
%定义边界条件=======================
delta=1;%用于存储两次计算的相对误差
%maxd=1;%N*N个相对误差中最大的一个
toi=0;%times of iteration迭代计算次数
while delta>tol
for m=2:N-1 n=2:N-1;
a(m,n)=u(m,n);
u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4;
delta=max(abs(u(m,n)-a(m,n))/u(m,n));%注意这种形式的意义。
end
toi=toi+1;
%maxd=max(delta(:));
end
X=1:N;Y=1:N;
u
mesh(X,Y,u(X,Y))
delta<tol
toi
toc
另外，在进行判断时对delta进行的比较也可这样编写：while delta>tol
delta=0;
for m=2:N-1 n=2:N-1;
a(m,n)=u(m,n);
u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4;
delta=max(abs(u(m,n)-a(m,n))/u(m,n),delta);
end
toi=toi+1;
%maxd=max(delta(:));
end
附：找矩阵中的最大值及其位置。
分情况（一维、二维或三维）而言：
i.一维阵，[a,b]=max(A)即可，a为最大值，b为位置;
ii.二维矩阵
a=max(A(:));
[x,y]=find(A==a);
iii.三维
a=max(A(:));
ind = find(a==max(a(:)));
[x,y,z] = ind2sub([m n d],ind);。