邊值問(wèn)題的研究.docx

西北農(nóng)林科技大學(xué)理學(xué)院應(yīng)用數(shù)學(xué)系微分方程數(shù)值解結(jié)課論文論文題目 邊值問(wèn)題的研究2016年 1 月 14 日一問(wèn)題重述對(duì)于下列邊值問(wèn)題：其中A為學(xué)號(hào)的倒數(shù)第2位，B為學(xué)號(hào)的倒數(shù)第1位。(1)差分：截?cái)嗾`差、穩(wěn)定性、收斂半徑、遞推（隱式）或方程組（顯式）(2)有限元：剛度矩陣、算法步驟及代碼二問(wèn)題分析題目明確指出使用差分方法和有限元解法。什么都不管先構(gòu)造一種差分格式，然后對(duì)求解區(qū)域做劃分將問(wèn)題離散化，從微分方程的定解問(wèn)題轉(zhuǎn)化為求線性代數(shù)方程的解，以便于能夠使用計(jì)算機(jī)進(jìn)行計(jì)算。在這里選用的是中心差分法，同時(shí)將邊界進(jìn)行處理，同時(shí)用Ritz有限元法和Galerkin法有限元法嘗試去得到結(jié)果，最后再去比較兩種解法所得到結(jié)果的精確性，分析相容性和截?cái)嗾`差等等。三解題過(guò)程1首先建立差分格式，考慮兩點(diǎn)的邊值問(wè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)離散化. 逼近邊值問(wèn)題（1.1）（1.2）的差分方程為：當(dāng)網(wǎng)格均勻，即時(shí)差分方程簡(jiǎn)化為這相當(dāng)于用一階中心差商，二階中心差商依次代替（1.1）的一階微商和二階微商的結(jié)果。這個(gè)方程就是中心差分格式。式（1.4）用方程組展開(kāi)：這是一個(gè)以為未知量的線性方程組。到此為止，中心差分格式展開(kāi)完畢，接下來(lái)處理方程(1.1)將方程在節(jié)點(diǎn)離散化，由泰勒公式展開(kāi)得：所以截?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代表邊界點(diǎn)，代表邊界點(diǎn)的取值。上式給出了差分方程的解的誤差估計(jì)，而且表明當(dāng)差分解收斂到原邊值問(wèn)題的解，收斂速度為。2接下來(lái)是有限元的解法從Ritz法出發(fā)，單元?jiǎng)偠染仃嚍椋喊匆?guī)則組裝成總剛度矩陣。令其中以及則有限元方程為從Galerkin有限元法出發(fā)，Galerkin有限元方程為：系數(shù)矩陣第j行只有三個(gè)非零元素，即這里第一行只有兩個(gè)非零元素：第n行只有兩個(gè)非零元素：和方程的右端項(xiàng)四求解過(guò)程其精確解為。算例中。(1)從Ritz法出發(fā)以將積分區(qū)間等分為10份為例，則步長(zhǎng)，記為。 為：以步長(zhǎng)取為h=1/10為例，從Ritz法出發(fā)的有限元法得到的數(shù)值解與精確解為Ritz數(shù)值解精確解001.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)邊值問(wèn)題很接近精確解。以步長(zhǎng)取為h=1/50為例，從Ritz法出發(fā)的有限元法得到的數(shù)值解與精確解圖像為最大誤差為0.002025步長(zhǎng)1/101/501/1001/500Ritz最大誤差0.20250.0441000.020250.004491分析：最大誤差為0.02025，Ritz有限元法求解兩點(diǎn)邊值問(wèn)題很接近精確解，且步長(zhǎng)越大，誤差越小。（2）從Galerkin法出發(fā)以將積分區(qū)間等分為10份為例，則步長(zhǎng)，記為。 為：Galerkin有限元法最大誤差:0.815000，圖像為： 以將積分區(qū)間等分為100份為例，圖像為：分析：最大誤差為0.080150，Galerkin有限元法求解兩點(diǎn)邊值問(wèn)題很接近精確解，且步長(zhǎng)越大，誤差越小。步長(zhǎng)1/101/501/1001/500Calerkin最大誤差0.801500.160060.0801500.016006最后收斂性和誤差分析：令和分別表示精確解和有限元解在剖分區(qū)間節(jié)點(diǎn)處的值，收斂性表示為記最大誤差為err,則問(wèn)題轉(zhuǎn)化為在方程中，已知h和err，求解M和k的擬合問(wèn)題。在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; %左邊界條件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); 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,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)=(subs(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,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=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.-,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_Galerkin);子程序：% Galerkin有限元法求解一維問(wèn)題%構(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); 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,0,1);fun4_Galerki