某二階偏微分方程計算方法對比分析

賈濤 牛華東 解景濤 盧瑞軍 蘇茂輝

摘要：CAE在汽車和發動機設計中發揮越來越重要的作用，然而CAE的本質是求解各種常微分方程組或偏微分方程組;常微分方程又稱動力系統方程，通常用來求解動力系統問題。偏微分方程則廣泛應用于電磁學、流體力學、結構力學等多個領域;對于偏微分方程通常很難求得其解析解，需要借助數值計算來獲取其方程的近似解。拉普拉斯方程（Laplace）是一種橢圓形二階偏微分方程，并且可以求取其解析解。本文通過解析法以及數值法對拉普拉斯方程求解，并對比不同求解方法的效率和精度;結論顯示解析法雖然精度較高，但是需要很大的計算量，并且大多數偏微分方程沒有解析解。因此，在汽車和發動機等工程應用中應該根據精度需求選擇最優的途徑求解偏微分方程問題。

關鍵詞：Laplace方程;偏微分方程;數值計算;CAE

1? 理論分析

拉普拉斯方程的數學表達式如下，

2? 問題描述

研究某矩形二維傳熱問題，已知傳熱問題的導熱微分方程式如下[2]：

3? 解析解

求解偏微分方程的核心在于邊界條件，不同的邊界條件甚至會出現完全不同的解;根據該問題描述，通過邊界條件將上述方程簡化：

由于公式（5）比較繁瑣，并且涉及到n（0～∞）求和計算;方程組規模小則計算量小，可以快速得到計算結果;方程組規模大則計算量很大，通常需要很久才能完成計算。本例中算法1取n=110，對5000個數據編程計算，則耗時近1小時;而采用算法2則只需要數秒時間。

算法1：

clc;

clear;

sum=zeros（100，50）;

syms x;

for i=1：100;

for j=1：50;

for k=1：110;

Dn=int（100*sin（k*pi*x/50），0，50）;

An=1/（25*sinh（2*pi*k））;

Cn=sinh（k*pi*（100-i）/50）;

Bn=sin（k*pi*j/50）;

Un=Dn*An*Cn*Bn;

sum（i，j）=sum（i，j）+Un;

end

end

end;

算法2：

clc;

clear;

sum=zeros（100，50）;

for i=1：100;

for j=1：50;

for k=1：110;

Dn=5000*（1-cos（k*pi））/（k*pi）;

An=1/（25*sinh（2*pi*k））;

Cn=sinh（k*pi*（100-i）/50）;

Bn=sin（k*pi*j/50）;

Un=Dn*An*Cn*Bn;

sum（i，j）=sum（i，j）+Un;

end

end;

end;

4? 數值解法

通過數值方法求解偏微分方程是常用的求解方法，如下通過有限元方法以及EXCEL迭代、軟件編程等不同方法進行數值計算，對比不同計算方法的精度和計算效率。

4.1 有限元方法

借助有限元軟件，建立平面模型并搭建穩態熱力學模型，生成5151個節點、10000個單元;設置邊界條件如上所述，計算得到的溫度場分布如圖2所示：該方法更偏向于工程的應用。在求解過程中，模型搭建過程花費較多時間;但是求解計算效率較高，普通電腦計算只需要20sCPU時間。

4.2 Excel迭代方法

設置步長dx=dy=1，采用五點差分法求解;如圖4所示，使用EXCEL建立數值計算表格，啟用迭代計算，設置誤差以及最大迭代次數進行求解;該方法計算效率較高，只需要1～2s CPU時間;得到計算結果如圖5[3]。

4.3 軟件編程計算

借助軟件編寫五點差分程序，使用高斯賽德爾迭代法，同樣可以計算出最終的結果;編程耗時較長，計算耗時60s左右的CPU時間[4]。

計算程序：

clear

clc

u=1：1：50;

v=1：1：50;

[x，y]=meshgrid（u，v）;%通過坐標軸點在平面畫網格

p=[x，y];

plot（p）;

spy（p）;

for i=2：50

for j=2：100

p（i，j）=（j-1）+98*（i-2）;

end;

end;

for i=1;

for j=1：100

p（i，j）=0;

end;

end;

for i=50;

for j=1：100

p（i，j）=0;

end;

end;

for j=1;

for i=1：50

p（i，j）=0;

end;

end;

for j=100;

for i=1：50

p（i，j）=0;

end;

end;

PP=delsq（p）;

load（'my.dat'）;

b=spconvert（my）;

N=length（b）;? ?%Gauss – Seidel迭代法

u=inv（PP）*b;

u=zeros（N，1）;%迭代初始值

D=diag（diag（PP））;

E=-tril（PP，-1）;%下三角

F=-triu（PP，1）;%上三角

B=inv（D-E）*F;g=inv（D-E）*b;

eps=0.000001;

for k=1：10000 %最大迭代次數

y=B*u+g;

if norm（u-y）<eps

break;

end

u=y;

end

s=zeros（50，100）;

for i=1：48

for j=1：98;

s（i+1，j+1）=u（98*（i-1）+j，1）;

end

end

for i=1：50

s（i，1）=100;

end

subplot（1，1，1）;

mesh（s）;

5? 誤差分析

對比計算結果第二列（x=1，u（x，y））數據并繪制圖表;如圖8所示，數值解和解析解的四條曲線重合在一起，說明這些數值大小相同，整體趨勢一致。圖9誤差分析顯示，EXCEL和編程計算誤差重合在一起，誤差在0～0.01%之間;有限元計算在中間部分誤差接近0，說明其計算誤差相對較小;三種數值計算方法誤差在工程可以接受的范圍內。

6? 總結

本文通過不同方法求解偏微分方程，對偏微分方程的求解思路有了深入的理解;各種解法中，效率最高的是使用EXCEL迭代求解，操作簡單并且計算速度快、精度較高。借助軟件編程的方法求解效率相對較低，并且需要對程序調試，精度方面和EXCEL求解差別不大。借助有限元求解則需要搭建有限元模型，并設置邊界條件，計算效率同樣很高，精度方面可以滿足工程需要。而解析法求解拉普拉斯方程雖然精度高，但是對算法要求較高。因此，對于偏微分方程工程方面應用，在滿足精度的前提下，解析解不一定是最優的方法;具體需要根據實際情況選擇相應的求解方法。

參考文獻：

[1]陸金甫，關治.偏微分方程數值解法[M].清華大學出版社，2004.

[2]楊世銘，陶文銓.傳熱學[M].四版.高等教育出版社，2006.

[3]賈新民，嚴文.有限差分法求解拉普拉斯方程[J].昌吉學院學報，2009.

[4]劉大衛，高明.正方形環域Laplace方程的簡明數值解法[J].沈陽師范大學學報，2006.