1、守而大修CENTRAL SOUTH UN IVERS ITY數值分析論文NUMERI CAL ANALYSIS (THESIS)題 目數值分析論文學生姓名指導教師張鴻雁學 院專業班級2016年12月目錄 TOC o 1-5 h z HYPERLINK l bookmark8 o Current Document 一二分法與牛頓迭代法1 HYPERLINK l bookmark30 o Current Document 二拉格朗日插值法6 HYPERLINK l bookmark38 o Current Document 三最小二乘法10 HYPERLINK l bookmark42 o Cur

2、rent Document 四復合辛普森求積公式15 HYPERLINK l bookmark57 o Current Document 五改進歐拉公式18 HYPERLINK l bookmark65 o Current Document 六列主元消去法21 HYPERLINK l bookmark71 o Current Document 七不動點迭代法26二分法與牛頓迭代法問題簡介在有氣體參加的恒容反應體系里，反應體系的總壓會隨著反應的進行而改變。通 過反應的平衡常數可以求出反應后氣體的總量，而氣體分子的分壓與其所占比率是成 正比的。這類問題中關鍵是要求出反應后各分子的分壓，各分子分壓之

3、和即為總壓，其實 質即為求方程的根。例如：在298K下，化學反應2OF2=Ch+2F2的平衡常數為0.410atm,如在298K下將OF? 通入容器，當t=0時為latm,問最后總壓是多少？計算精度為IO%二數學模型假設是理想氣體，可知2OF2=O2+2F2設氧的分壓為p,平衡時有l-2p p p平衡時有-2=0.410(l-2p)2整理得4P3-1.640P2+ 1.64p-0.410=0函數關系式為f (p) = 4p3 - 1.640p2 + 1.64p - 0.410 = 0三算法選擇與算法過程由計算得f(0.2)=-0.1156, f(0.3尸0.0424因此，有根區間為0.2, 0

4、.3用求單根的二分法計算時：可令a=0.2, b=0.3；最后計算得氧氣分壓：p=0.274609atm,總壓為：P=3p+( 1 -2p)= 1.274609atm用牛頓迭代法計算時:令f(x) = 4x3 - 1.640 x2 + 1.64x - 0.410f(x) = 12x2 -3.280 x + 1.64迭代公式為：x+ = x-黑令x = 0.2進行迭代，最后計算得氧氣分壓為：p=0.274901,則總壓為:P=3p+(l-2p)= 1.274901二分法的計算算法:#include #includeusing namespace std;double a=0.2,b=0.3,e=

5、0.001,n;double f(float c)(n=4*c*c*c-1.6400000*c*c+1.64*c-0.410;retiun n;void mainO(double c,m,p,a=0.2,b=0.3;fbr(p=a;b-a=0.001;)(c=(a+b)/2;m=f(c);if(m0)b=c;elsep=c;break;p=(a+b)/2;coutn p=upendl;coutHn a=naendl;coutMn b=nbendl;牛頓迭代法的算法：#include#include#includeusing namespace std;double f(double x)(do

6、uble m,n;m=4*x*x*x-1.640*x*x+l .64*x-0.410;n=12*x*x-3.280*x+1.64;retxuii x-m/n;void newton(double x,double d)(double a=x;double b=f(a);hit k; /記錄循環的次數fbr(k=l ;fabs(a-b)d;k+)(a=b;b=f(a);if(k100)(coutvu迭代失敗，該函數可能不收斂!yvendl;coutn a=Haendl;coutnn b=ubendl;coutnn k=nkendl;retiun;int mainOCOUtVV”請輸入初始值XO和

7、要求得結果的精度：”;double x,d;cinxd;newton(x,d);retxuii 0;數值實駛過程二分法通過VC6.0程序運行的結果如下圖所示:牛頓迭代法通過VC6.0的計算結果如下:五相關數值分析和實際應用分析由二分法程序輸出結果可得:a=0.274219,b=0.275。所以實驗誤差： |x*-Xk|W (b-a) /2=0.00040.001滿足實驗要求。由牛頓迭代法程序輸出結果為a=0.274918, b=0.274901o實驗誤差:|x*-Xk|a-b|=0.0000170.001,滿足實驗要求。二拉格朗日插值法-問題簡介在化學實驗中，通常測得的是一批離散數據，需要從這

8、批有限的測量數據中得出 一個函數關系，進而來求解任意的函數值，如：實驗測得某物質在20c下，其粘度(Pa-s)與水溶液濃度c (重量)的關系如下 表。試用拉格朗日5次插值計算粘度在2.0X10-3和5.5X10-3時所對應的濃度。要求精 確到小數點后四位數。粘度（Pa-s）c （重量）粘度（Pa-s）c （重量）1.005 XI O-303.652XW3401.776X10-3204.621X10-3452.480 XI0-3305.921 XW350二數學模型在一定溫度下，物質的水溶液粘度與其濃度是呈正相關的，表中給出6個點的函數值，所以可以進行拉格朗日5次插值計算，其計算公式為:k=0其中

9、X即為粘度，k(x)即為所求值溶液濃度。三算法選擇與算法過程計算算法：#includeusing namespace std;const int N=1000;float Lagi,ange（float aiTX,float arrY口,int n,float x）float yResult=0.0;float LalueN;int k,m;float templ,temp2;fbr(k=0;kn;k+)(temp 1=1.0;temp2=1.0;fbr(m=0 ;mn;m+)(if(m=k)(continue;)tempi *= (x-anXm);temp2 *= (aiTXk-aiTXm)

10、;Laluek=temp l/temp2;fbr(int i=0;in;i+)(yResult += aiTYi*Laluei;retiun yResult;)int maiii()(float aiTXN,an-YN;int num;coutvv”輸入插值節點的個數(小于“vN個):H；cinnum;coutvv”n一接下來輸入這些插值節點(先輸入X再輸入對應的Y)-n”;fbr(int i=O;inum;i+)coutw”第Yvi+iw”個節點的X值：H;cinaiTXi;coutvv第“vvi+lvv”個節點的 Y 值:； cinarrYi;float X;coutvv”n-請輸入第一個

11、待求解的插值節點的X值-cinX;float Resl = Lagraiige(anX,aiTY；num,X);coutvv”n一第一個點插值結果為：MRes lendl;coutvv”n請輸入第二個待求解的插值節點的X值一n”;cinX;float Res2 = Lagrange(arrX,anYnum,X);coutvv”n一第二個點插值結果為：,Res2endl;retiun 0;四數值實駛過程Ax YXYXYXYXYXY 耕點占或點點點點點點點點點 局節節節節節節節節節節節節 書 112233445566/ y I y i y i 1/ y I Jz y / i y i Jzxyxyx

12、v-xvxvxy一第第第第第第第第第第第第五相關數值分析和實際應用分析在實驗過程中，拉格朗日高次差值對數據的誤差很敏感，由于舍入誤差，在計算 中這種誤差很可能被放大，因此在實際工作中，節點數一般不超過七個。最小二乘法-問題簡介用紫外吸收分光光度計測定牛血清蛋白的紫外吸收標準曲線時，先配制LOmg/mL 的牛血清蛋白溶液，然后分別稀釋為O.125mg/mL, O.25mg/mL, O.375mg/mL,O.5mg/mL, 0.625mg/mL, O.75mg/mL, O.875mg/niL的溶液，分別測定不同濃度溶液在278mn處 的紫外吸收，實驗數據記錄如下表，試利用最小二乘法計算出濃度c與吸

13、光度A之間 的關系。cAcA0.1250.0620.6250.3220.250.1240.750.3830.3750.1900.8750.4520.50.25510.512二數學模型物質的紫外吸收吸光度與它的濃度滿足朗伯比爾定律，即：A= Kbc式中A為吸光度，K為摩爾吸收系數，它與吸收物質的性質及入射光的波長人有 關，c為吸光物質的濃度，b為吸收層厚度。當保證入射光波長與吸收層厚度不變時， 吸光度A與濃度c是呈線性相關的。三算法選擇與算法過程選用最小二乘法做線性擬合。計算算法：#inc lude#includedouble f(double a, int b)int i;double k=

14、1;for (i=0; ib; i+)retxuii k;int main(void)(int i, j, n, m;double a22, d2=0, e2, t;printf(請輸入Xi和Yi的個數：n)；scanf(n%dM, &n);double *b = (double *) malloc(sizeof(double)*n);double *c = (double *) malloc(sizeof(double)*n);printf(”請輸入 Xi： nH);for (i=0; in; i+)scaiif(,%lf, &bi);printf(”請輸入 Yi： nH);fbr (j=0

15、; jn; j+)scaiif(,%lf, &cj);fbr (i=0; i2; i+)(for (j=0; j2;j+)(aiD = O；for (m=0; mn; m+)aiU += f(bm,i*2+j);all = alO;alO = aOl;for (i=0; in; i+)(d0 += ci;dl+=bi*ci;fbr(i=O; i2; +i)(t = aOi; aOi = ali;ali = t;t = dO;dO = dl;dl=t;all-=aOl*alO/aOO;dl -=dO*alO/aOO;alO = O;if(all = 0 & dl !=0)printf(“無解。)

16、；else(el = dl/all;eO = (dO-el*aOl)/aOO;printf(方程的系數為:nn);for (j=0; j2;j+)printf(A%d = %10.61fiiM,j+l, ej);printf(方程為:n)printf(MP(x) = %10.61f + %10.61fxir, e0, el);fiee(b);fiee(c);四數值實驗過程五相關數值分析和實際應用分析由表中數據畫出的A-c的圖像如下圖所示:在分析化學工作中，物質的濃度與溶液的吸光度、電位或其它物理性質存在線性 關系，因此在測量中經常要計算標準工作曲線法，而最小二乘法是在求標準曲線時的 一種方便有

17、效的方法。四復合辛普森求積公式-問題簡介用脈沖法測定某反應器的停留時間分布，由實驗數據計算壽命分布密度與時間t 的關系如下表t10152025303540E(t)X10-301.929.020.9531.1032.2030.60求分布密度積分：優5 E(t)dt二數學模型根據所給數據，小區間分為4等分，步長h=10, n=4,用復合辛普森求積公式計 算。 TOC o 1-5 h z n-1n-1Sn = J f (q) + f + 4 W f (，+J + 2 W f(Xk)k=0k=lf(a尸f(10尸0 f(b)=f(40)=30.60 X IO3dt-O + 30.6 +4(1.92 +

18、 20.95 + 32.2) + 2(9.0 + 31.1) X IO-3 = 0.55181U6三算法過程由于本題中的函數值都為散點，令m=2k,此時n=7, h=5則：n-1n-1Sn=gr(a)+r(0) + 4 W f(Xm)+2 W f(%m)m=l,3,.m=2,4,.計算算法：#includeusing namespace std;double Xiiipusen(double a,double b,int n)intdouble y7=0,1.92,9.0,20.95,31.10,32.20,30.60;double h=(b-a)/n;double s=y0+y6;for(k

19、=0;kn/2;k+)(m=2*k+l;s+=4*ym;fbr(k=l ;kii/2;k-H-)(t=2*k;s+=2*yt;retxuii s*0.001*lV3;int mainO(double a=10;double b=40;int n=6;coutM由復化辛普森公式求的結果：vvXinpusen(a,b,n)wendl： retxuii 0;四數值實驗過程*D:C-h+MSDev98MyProjectsxpsDebugxps.exe,由復化辛苦森公式求的事 0.5518Press any key to continue五相關數值分析和實際應用分析本實驗通過10到40分鐘的反應壽命分布

20、密度數據來計算反應器的停留時間分布, 利用復合辛普森求積法對10到40分鐘內的壽命分布密度進行積分，即得10到40分 鐘內反應停留時間占總反應時間的比例。實驗計算結果為0.5518,即10到40分鐘內 的反應停留時間占總反應時間的比率為55.18%。五改進歐拉公式一問題簡介已知在管式反應器中進行液相反應A-R+S,反應為吸熱反應，反應管外油浴溫 度為340C,假定已知管內溫度與轉化率的關系為：dt=-65.0 - 15.58(1 - %)/(伙1 -4)速度常數為：k=1.17X 1017exp(-22400/T)tc=反應器外壁溫度。如反應器入口溫度為340,反應器出口轉化率為90%,試用改

21、 進歐拉公式求反應器出口溫度。出口 =轉化率、=90%，t=? 甘入口 ：轉化率t=340 二數學模型由題意知為常微分方程的初值問題，假定=340(管外油浴溫度)，方程形式為:dt-=-65.0 - 15.58(t - 340)/(k(l - xA) k=1.17X 1017exp(-22400/(t+273.2)Xa=0，to=34O使用改進歐拉法計算，設自變量xa的步長h=0.05,經過18步，可求出xa=0.9時的出口溫度。其數學模型為:2xnyp = % + h(% -)Jn2xn+1yc =yn +h(# )yp1yn+i = 5(4 + %)三算法選擇與算法過程使用改進歐拉法計算，

22、設自變量xa的步長h=0.05,經過18步，可求出xa=0.9時的出口溫度。程序分為兩部分：(a)定義方程右邊函數；(b)用改進歐拉法計算。計算算法：#include#include#define N 18void ModEuler(float (*f 1 )(float,float),float xO,float yO.float xn,int n)(int i;float yp,yc,x=xO,y=yO,h=(xn-xO)/n;cout,x0=,x,ty0,yendl;fbr(i=l;i x3,查得Na2c03、NaHCO3和NaOH的分子量分別為100、80、40。由題意可列出如下線性方

23、程組:+ %2 + X3 = 200 x3房 + :=0.1 X 8100 402%1 %2%3奈+關+2=0.1 X 30100 80 401可令 A= 1/1002/100101/4011/401/40200 b = 0.83則:Ax=b三算法選擇與算法過程選用列主元消去法解上述線性方程組。1A = 0.010.02100.012510.0250.025計算算法:#include#include#includeusing namespace std;#defiiie maxn 50int n;double a maxn maxn;double bmaxn;/double xmaxn;/vo

24、id read()(coutvv”請輸入系數矩陣規模n:= ； cinn;coutw”;請輸入系數矩陣：nu;for(int i=l;i=n;i+)for(intj=l;j=n;j+) cinaij;)coutw請輸入b矩陣：n”；for(i=l ;i=ii;i-H-) cinbi;void PrintQ(cout|nu;coutv”|結果為：n”；for(int i=l;i=n;i+)(cout,x,i,=,xiendl;cout|nnH;void LieZhuXiaoYtian()for(intk=l;kn;k+)double ab_max=-l;mt max_ik;for(int i=k

25、;iab_max)(ab_max=abs(a ik); max_ik=i;)if(max_ik!=k)(double temp;fbr(intj=l;j=ii;j-H-)(temp=a max_ik j ;amax_ikj=akj;akj=temp;temp=bmax_ik;bmax_ik=bk;bk=temp;)fbr(i=k+1 ;i=ii;i+)aik/=akk;fbr(int j =k+1 ;j =ii;j +)aiU-=aik*akU;bi-=aik*bk;)if(kO;i)(xi=bi；for(int j=i+l;j=n;j+)(xi-=aiU*xU； xi/=aii；Print(

26、);)int main()while(l)read();LieZhuXiao Yiian();四數值實駛過程五 相關數值分析和實際應用分析高斯消去法既容易以理解乂很簡單，但經常得不到正確結果，原因是消元過程是 按方程組的自然順序進行的。如果消去元數值太小，會產生較大的舍入誤差，得出完 全錯誤的結果。如果消去元是零會導致消元過程無法進行。為了解決這些問題，就要 選用列主元消去法。在根據光度確定多組分溶液中各組分濃度時也常常用到列主元消 去法。七不動點迭代法一問題簡介在9.33atm, 302K時，容器中充2moi氨氣，試求該容器的容積，要求精度為IO, （氨氣的范德華常數:a=4.17atm L2/mol2,b=0.0371 atm L2/mol2 ）二數學模型范德華方程:改寫為迭代公式并代入數據:nRTV =y