數(shù)值分析上機(jī)C++版

1、守而大修CENTRAL SOUTH UN IVERS ITY數(shù)值分析論文NUMERI CAL ANALYSIS (THESIS)題 目數(shù)值分析論文學(xué)生姓名指導(dǎo)教師張鴻雁學(xué) 院專業(yè)班級(jí)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 四復(fù)合辛普森求積公式15 HYPERLINK l bookmark57 o Current Document 五改進(jìn)歐拉公式18 HYPERLINK l bookmark65 o Current Document 六列主元消去法21 HYPERLINK l bookmark71 o Current Document 七不動(dòng)點(diǎn)迭代法26二分法與牛頓迭代法問題簡介在有氣體參加的恒容反應(yīng)體系里，反應(yīng)體系的總壓會(huì)隨著反應(yīng)的進(jìn)行而改變。通 過反應(yīng)的平衡常數(shù)可以求出反應(yīng)后氣體的總量，而氣體分子的分壓與其所占比率是成 正比的。這類問題中關(guān)鍵是要求出反應(yīng)后各分子的分壓，各分子分壓之

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

4、.3用求單根的二分法計(jì)算時(shí)：可令a=0.2, b=0.3；最后計(jì)算得氧氣分壓：p=0.274609atm,總壓為：P=3p+( 1 -2p)= 1.274609atm用牛頓迭代法計(jì)算時(shí):令f(x) = 4x3 - 1.640 x2 + 1.64x - 0.410f(x) = 12x2 -3.280 x + 1.64迭代公式為：x+ = x-黑令x = 0.2進(jìn)行迭代，最后計(jì)算得氧氣分壓為：p=0.274901,則總壓為:P=3p+(l-2p)= 1.274901二分法的計(jì)算算法:#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; /記錄循環(huán)的次數(shù)fbr(k=l ;fabs(a-b)d;k+)(a=b;b=f(a);if(k100)(coutvu迭代失敗，該函數(shù)可能不收斂!yvendl;coutn a=Haendl;coutnn b=ubendl;coutnn k=nkendl;retiun;int mainOCOUtVV”請(qǐng)輸入初始值XO和

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

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

9、X即為粘度，k(x)即為所求值溶液濃度。三算法選擇與算法過程計(jì)算算法：#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”輸入插值節(jié)點(diǎn)的個(gè)數(shù)(小于“vN個(gè)):H；cinnum;coutvv”n一接下來輸入這些插值節(jié)點(diǎn)(先輸入X再輸入對(duì)應(yīng)的Y)-n”;fbr(int i=O;inum;i+)coutw”第Yvi+iw”個(gè)節(jié)點(diǎn)的X值：H;cinaiTXi;coutvv第“vvi+lvv”個(gè)節(jié)點(diǎn)的 Y 值:； cinarrYi;float X;coutvv”n-請(qǐng)輸入第一個(gè)

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

12、v-xvxvxy一第第第第第第第第第第第第五相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析在實(shí)驗(yàn)過程中，拉格朗日高次差值對(duì)數(shù)據(jù)的誤差很敏感，由于舍入誤差，在計(jì)算 中這種誤差很可能被放大，因此在實(shí)際工作中，節(jié)點(diǎn)數(shù)一般不超過七個(gè)。最小二乘法-問題簡介用紫外吸收分光光度計(jì)測定牛血清蛋白的紫外吸收標(biāo)準(zhǔn)曲線時(shí)，先配制LOmg/mL 的牛血清蛋白溶液，然后分別稀釋為O.125mg/mL, O.25mg/mL, O.375mg/mL,O.5mg/mL, 0.625mg/mL, O.75mg/mL, O.875mg/niL的溶液，分別測定不同濃度溶液在278mn處 的紫外吸收，實(shí)驗(yàn)數(shù)據(jù)記錄如下表，試?yán)米钚《朔ㄓ?jì)算出濃度c與吸

13、光度A之間 的關(guān)系。cAcA0.1250.0620.6250.3220.250.1240.750.3830.3750.1900.8750.4520.50.25510.512二數(shù)學(xué)模型物質(zhì)的紫外吸收吸光度與它的濃度滿足朗伯比爾定律，即：A= Kbc式中A為吸光度，K為摩爾吸收系數(shù)，它與吸收物質(zhì)的性質(zhì)及入射光的波長人有 關(guān)，c為吸光物質(zhì)的濃度，b為吸收層厚度。當(dāng)保證入射光波長與吸收層厚度不變時(shí)， 吸光度A與濃度c是呈線性相關(guān)的。三算法選擇與算法過程選用最小二乘法做線性擬合。計(jì)算算法：#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(請(qǐng)輸入Xi和Yi的個(gè)數(shù)：n)；scanf(n%dM, &n);double *b = (double *) malloc(sizeof(double)*n);double *c = (double *) malloc(sizeof(double)*n);printf(”請(qǐng)輸入 Xi： nH);for (i=0; in; i+)scaiif(,%lf, &bi);printf(”請(qǐng)輸入 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(方程的系數(shù)為: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);四數(shù)值實(shí)驗(yàn)過程五相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析由表中數(shù)據(jù)畫出的A-c的圖像如下圖所示:在分析化學(xué)工作中，物質(zhì)的濃度與溶液的吸光度、電位或其它物理性質(zhì)存在線性 關(guān)系，因此在測量中經(jīng)常要計(jì)算標(biāo)準(zhǔn)工作曲線法，而最小二乘法是在求標(biāo)準(zhǔn)曲線時(shí)的 一種方便有

17、效的方法。四復(fù)合辛普森求積公式-問題簡介用脈沖法測定某反應(yīng)器的停留時(shí)間分布，由實(shí)驗(yàn)數(shù)據(jù)計(jì)算壽命分布密度與時(shí)間t 的關(guān)系如下表t10152025303540E(t)X10-301.929.020.9531.1032.2030.60求分布密度積分：優(yōu)5 E(t)dt二數(shù)學(xué)模型根據(jù)所給數(shù)據(jù)，小區(qū)間分為4等分，步長h=10, n=4,用復(fù)合辛普森求積公式計(jì) 算。 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三算法過程由于本題中的函數(shù)值都為散點(diǎn)，令m=2k,此時(shí)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,.計(jì)算算法：#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由復(fù)化辛普森公式求的結(jié)果：vvXinpusen(a,b,n)wendl： retxuii 0;四數(shù)值實(shí)驗(yàn)過程*D:C-h+MSDev98MyProjectsxpsDebugxps.exe,由復(fù)化辛苦森公式求的事 0.5518Press any key to continue五相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析本實(shí)驗(yàn)通過10到40分鐘的反應(yīng)壽命分布

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

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

22、設(shè)自變量xa的步長h=0.05,經(jīng)過18步，可求出xa=0.9時(shí)的出口溫度。程序分為兩部分：(a)定義方程右邊函數(shù)；(b)用改進(jìn)歐拉法計(jì)算。計(jì)算算法：#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奈+關(guān)+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計(jì)算算法:#include#include#includeusing namespace std;#defiiie maxn 50int n;double a maxn maxn;double bmaxn;/double xmaxn;/vo

24、id read()(coutvv”請(qǐng)輸入系數(shù)矩陣規(guī)模n:= ； cinn;coutw”;請(qǐng)輸入系數(shù)矩陣：nu;for(int i=l;i=n;i+)for(intj=l;j=n;j+) cinaij;)coutw請(qǐng)輸入b矩陣：n”；for(i=l ;i=ii;i-H-) cinbi;void PrintQ(cout|nu;coutv”|結(jié)果為：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();四數(shù)值實(shí)駛過程五 相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析高斯消去法既容易以理解乂很簡單，但經(jīng)常得不到正確結(jié)果，原因是消元過程是 按方程組的自然順序進(jìn)行的。如果消去元數(shù)值太小，會(huì)產(chǎn)生較大的舍入誤差，得出完 全錯(cuò)誤的結(jié)果。如果消去元是零會(huì)導(dǎo)致消元過程無法進(jìn)行。為了解決這些問題，就要 選用列主元消去法。在根據(jù)光度確定多組分溶液中各組分濃度時(shí)也常常用到列主元消 去法。七不動(dòng)點(diǎn)迭代法一問題簡介在9.33atm, 302K時(shí)，容器中充2moi氨氣，試求該容器的容積，要求精度為IO, （氨氣的范德華常數(shù):a=4.17atm L2/mol2,b=0.0371 atm L2/mol2 ）二數(shù)學(xué)模型范德華方程:改寫為迭代公式并代入數(shù)據(jù):nRTV =y