邊值問題的研究

1、西北農(nóng)林科技大學(xué)理學(xué)院應(yīng)用數(shù)學(xué)系微分方程數(shù)值解結(jié)課論文論文題目 邊值問題的研究2016年 1 月 14 日一·問題重述對(duì)于下列邊值問題：其中A為學(xué)號(hào)的倒數(shù)第2位，B為學(xué)號(hào)的倒數(shù)第1位。(1)差分：截?cái)嗾`差、穩(wěn)定性、收斂半徑、遞推（隱式）或方程組（顯式）(2)有限元：剛度矩陣、算法步驟及代碼二·問題分析題目明確指出使用差分方法和有限元解法。什么都不管先構(gòu)造一種差分格式，然后對(duì)求解區(qū)域做劃分將問題離散化，從微分方程的定解問題轉(zhuǎn)化為求線性代數(shù)方程的解，以便于能夠使用計(jì)算機(jī)進(jìn)行計(jì)算。在這里選用的是中心差分法，同時(shí)將邊界進(jìn)行處理，同時(shí)用Ritz有限元法和Galerkin法有限元法嘗試

2、去得到結(jié)果，最后再去比較兩種解法所得到結(jié)果的精確性，分析相容性和截?cái)嗾`差等等。三·解題過程1·首先建立差分格式，考慮兩點(diǎn)的邊值問題，由題目知道建立中心差分格式如下對(duì)求解區(qū)間做網(wǎng)格劃分，在a到b之間取N+1個(gè)節(jié)點(diǎn)，定義為xi（i取1到N）即將區(qū)間I=a,b分為N個(gè)小區(qū)間由此得到區(qū)間的一個(gè)網(wǎng)格剖分。記。用表示網(wǎng)格內(nèi)點(diǎn)，的集合，表示內(nèi)點(diǎn)和界點(diǎn)的集合。取相鄰節(jié)點(diǎn)的中點(diǎn)，稱為半整數(shù)點(diǎn)。由節(jié)點(diǎn)又構(gòu)成的一個(gè)對(duì)偶剖分。用差商代替微商，將方程（1.1）在內(nèi)點(diǎn)離散化. 逼近邊值問題（1.1）（1.2）的差分方程為：當(dāng)網(wǎng)格均勻，即時(shí)差分方程簡化為這相當(dāng)于用一階中心差商，二階中心差商依次代替（1.

3、1）的一階微商和二階微商的結(jié)果。這個(gè)方程就是中心差分格式。式（1.4）用方程組展開：這是一個(gè)以為未知量的線性方程組。到此為止，中心差分格式展開完畢，接下來處理方程(1.1)將方程在節(jié)點(diǎn)離散化，由泰勒公式展開得：所以截?cái)嗾`差為下一步是分析差分格式的穩(wěn)定性差分格式的截?cái)嗾`差：，而邊界條件的截?cái)嗾`差為收斂性和穩(wěn)定性是從不同角度討論差分法的精確情況，穩(wěn)定性主要是討論初值的誤差和計(jì)算中的舍入誤差對(duì)計(jì)算結(jié)果的影響，收斂性則主要討論推算公式引入的截?cái)嗾`差對(duì)計(jì)算結(jié)果的影響.使用既收斂有穩(wěn)定的差分格式才有比較可靠的計(jì)算結(jié)果，這也是討論收斂性和穩(wěn)定性的重要意義.截?cái)嗾`差：，即。差分方程組的解滿足：其中a、b代表邊

4、界點(diǎn)，代表邊界點(diǎn)的取值。上式給出了差分方程的解的誤差估計(jì)，而且表明當(dāng)差分解收斂到原邊值問題的解，收斂速度為。2·接下來是有限元的解法從Ritz法出發(fā)，單元?jiǎng)偠染仃嚍椋喊匆?guī)則組裝成總剛度矩陣。令其中以及則有限元方程為從Galerkin有限元法出發(fā)，Galerkin有限元方程為：系數(shù)矩陣第j行只有三個(gè)非零元素，即這里第一行只有兩個(gè)非零元素：第n行只有兩個(gè)非零元素：和方程的右端項(xiàng)四·求解過程其精確解為。算例中。(1)從Ritz法出發(fā)以將積分區(qū)間等分為10份為例，則步長，記為。 為：以步長取為h=1/10為例，從Ritz法出發(fā)的有限元法得到的數(shù)值解與精確解為Ritz數(shù)值解精確解00

5、1.15401.11352.22452.14803.20253.09454.07903.94404.84504.68755.49155.31606.00955.82056.39006.19206.62406.42156.702506.5000圖像為分析：最大誤差為0.202500，Ritz有限元法求解兩點(diǎn)邊值問題很接近精確解。以步長取為h=1/50為例，從Ritz法出發(fā)的有限元法得到的數(shù)值解與精確解圖像為最大誤差為0.002025步長1/101/501/1001/500Ritz最大誤差0.20250.0441000.020250.004491分析：最大誤差為0.02025，Ritz有限元法求解

6、兩點(diǎn)邊值問題很接近精確解，且步長越大，誤差越小。（2）從Galerkin法出發(fā)以將積分區(qū)間等分為10份為例，則步長，記為。 為：Galerkin有限元法最大誤差:0.815000，圖像為： 以將積分區(qū)間等分為100份為例，圖像為：分析：最大誤差為0.080150，Galerkin有限元法求解兩點(diǎn)邊值問題很接近精確解，且步長越大，誤差越小。步長1/101/501/1001/500Calerkin最大誤差0.801500.160060.0801500.016006最后收斂性和誤差分析：令和分別表示精確解和有限元解在剖分區(qū)間節(jié)點(diǎn)處的值，收斂性表示為記最大誤差為err,則問題轉(zhuǎn)化為在方程中，已知h和e

7、rr，求解M和k的擬合問題。在Matlab中擬合采用最小二乘法實(shí)現(xiàn)。對(duì)和進(jìn)行最小二乘冪函數(shù)擬合，求得從Ritz法出發(fā)的誤差階為k=0.9612,M=0.4115.對(duì)和進(jìn)行最小二乘冪函數(shù)擬合，求得從Galerkin法出發(fā)的誤差階為k=1.004,M=3.061.五·操作代碼主程序：function Ritz(a,b,N)% -D2u=a*x+b，u(0)=0，Du(1)=0%a=9;b=7; %a為學(xué)號(hào)倒數(shù)第二位，b為學(xué)號(hào)倒數(shù)第一位，%N為剖分份分?jǐn)?shù)% 調(diào)用格式：Ritz(9,7,10); %N為剖分份分?jǐn)?shù)syms x;ua=0; %區(qū)間界點(diǎn)ub=1; %區(qū)間界點(diǎn)u_a=0; %左邊界

8、條件du_b=0; %右邊界條件p=1;q=0;f=a*x+b; %f=9x+7h=1/N;X=0:h:1; K=Stiff_matrix(p,q,h,N,X); %得到總剛度矩陣KB=vector(p,q,X,h,N,f); %得到BU = KB; %差分解 %處理右邊值條件u_b = (2*h*du_b-U(end-1)+4*U(end)/3;U=0;U;u_b; %精確解y0 = dsolve('-D2y=9*X+7','y(0)=0','Dy(1)=0','X');precise_value=double(subs(y0)

9、; deta=U-precise_value' deta_max=max(abs(deta); %最大誤差fprintf('最大的誤差是%fn',deta_max)% 差分解與精確解對(duì)比表figure;subplot(1,2,1);plot(X,U,'b*',X,precise_value,'r-');xlabel('x');ylabel('u');title('差分解與精確解對(duì)比圖');legend('差分解','精確解','Location'

10、;,'best');% 誤差圖subplot(1,2,2);plot(X,deta);xlabel('x');ylabel('u');end子程序：得到剛度矩陣K：function K=Stiff_matrix(p,q,h,N,X) syms x;K=zeros(N-1,N-1); diag_0=zeros(N-1,1); %確定K的對(duì)角元diag_1=zeros(N-2,1); %確定K的偏離對(duì)角的上對(duì)角元diag_2=zeros(N-2,1); %確定K的偏離對(duì)角的下對(duì)角元 % 獲取對(duì)角元素for j=1:N-1 diag_0(j)=(sub

11、s(p,x,(X(j)+X(j+1)/2)+(subs(p,x,(X(j)+X(j+1)/2+h)+(subs(q,x,X(j+1)*(h2);end% 獲取A的第三條對(duì)角for j=1:N-2 diag_2(j)=-(subs(p,x,(X(j+1)+X(j+2)/2)-(subs(0,x,X(j+2)*h)/2;end%獲取A的第二條對(duì)角for j=1:N-2 diag_1(j)=-(subs(p,x,(X(j+1)+X(j+2)/2)+(subs(0,x,X(j+1)*h)/2;endK=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1); K(N-1

12、,N-1)=2;K(N-1,N-2)=-2;得到B：function B=vector(p,q,x,h,N,f)B=zeros(N-1,1); syms x;X=0:h:1;for j=1:N-1; B(j)=(h2)*(subs(f,x,X(j+1);end%系數(shù)矩陣B(1)=B(1)+0*(subs(p,x,(X(2)+X(1)/2)+(subs(0,x,X(2)*h/2);B(N-1)=3*B(N-1)+2*0*(subs(p,x,(X(N)+X(N+1)/2)-(subs(0,x,X(N)*h/2)主程序：p=1;q=0;a = 9;b = 7;N = 100;%剖分份數(shù)h=1/N;x

13、=0:h:1;A_Galerkin=matirix(a,b,p,q,h,N);values_f_Galerkin=vector1(a,b,x,h,N);U_Galerkin=A_Galerkinvalues_f_Galerkin; y0 = dsolve('-D2y=a*x+b','y(0)=0','Dy(1)=0','x');precise_value=double(subs(y0);% Galerkin有限元法解與精確解對(duì)比圖figure;%subplot(1,2,1);plot(x,0;U_Galerkin,'b.-

14、',x,precise_value,'r.-');xlabel('x');ylabel('u');title('Galerkin有限元法解與精確解對(duì)比圖');legend('Galerkin數(shù)值解','精確解','Location','best');err_Galerkin=max(abs(0;U_Galerkin-precise_value');fprintf(sprintf('Galerkin有限元法最大誤差:%fn',err_Ga

15、lerkin);子程序：% Galerkin有限元法求解一維問題%構(gòu)造系數(shù)矩陣function A_Galerkin=matirix(a,b,p,q,h,N),A_Galerkin=zeros(N);for i=3:N fun1_Galerkin=(ks)-p./h+h.*q.*ks.*(1-ks); fun2_Galerkin=(ks)p./h+h.*q.*(ks.2)+p./h+h.*q.*(1-ks).2); fun3_Galerkin=(ks)-p./h+h.*q.*ks.*(1-ks); A_Galerkin(i-1,i-2)=integral(fun1_Galerkin,0,1);

16、 A_Galerkin(i-1,i-1)=integral(fun2_Galerkin,0,1); A_Galerkin(i-1,i) =integral(fun3_Galerkin,0,1);endfun2_Galerkin=(ks)p./h+h.*q.*(ks.2)+p./h+h.*q.*(1-ks).2);A_Galerkin(1,1)=integral(fun2_Galerkin,0,1);fun3_Galerkin=(ks)-p./h+h.*q.*ks.*(1-ks);A_Galerkin(1,2)=integral(fun3_Galerkin,0,1);fun3_Galerkin=(ks)-p./h+h.*q.*ks.*(1-ks);A_Galerkin(N,N-1)=integral(fun3_Galerkin