稳定温度场的拉普拉斯方程

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

toi

toc