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

1、 數(shù)值分析論文NUMERICAL ANALYSIS(THESIS)題 目數(shù)值分析論文學(xué)生姓名指導(dǎo)教師張鴻雁學(xué) 院專(zhuān)業(yè)班級(jí)2016年12月目錄一 二分法與牛頓迭代法1二 拉格朗日插值法 6三 最小二乘法10四 復(fù)合辛普森求積公式15五 改進(jìn)歐拉公式18六 列主元消去法21七 不動(dòng)點(diǎn)迭代法26一 二分法與牛頓迭代法一 問(wèn)題簡(jiǎn)介在有氣體參加的恒容反應(yīng)體系里，反應(yīng)體系的總壓會(huì)隨著反應(yīng)的進(jìn)行而改變。通過(guò)反應(yīng)的平衡常數(shù)可以求出反應(yīng)后氣體的總量，而氣體分子的分壓與其所占比率是成正比的。這類(lèi)問(wèn)題中關(guān)鍵是要求出反應(yīng)后各分子的分壓，各分子分壓之和即為總壓，其實(shí)質(zhì)即為求方程的根。例如：在298K下，化學(xué)反應(yīng)2OF2

2、=O2+2F2的平衡常數(shù)為0.410atm，如在298K下將OF2通入容器，當(dāng)t=0時(shí)為1atm，問(wèn)最后總壓是多少？計(jì)算精度為10-3。二 數(shù)學(xué)模型假設(shè)是理想氣體，可知 2OF2=O2+2F2設(shè)氧的分壓為p，平衡時(shí)有 1-2p p p平衡時(shí)有 4p3(1-2p)2=0.410整理得 4p3-1.640p2+1.64p-0.410=0函數(shù)關(guān)系式為 fp=4p3-1.640p2+1.64p-0.410=0三 算法選擇與算法過(guò)程由計(jì)算得 f(0.2)=0.1156，f(0.3)=0.0424因此，有根區(qū)間為0.2，0.3用求單根的二分法計(jì)算時(shí)：可令a=0.2，b=0.3；最后計(jì)算得氧氣分壓：p=0.

3、274609atm，總壓為：P=3p+(1-2p)=1.274609atm用牛頓迭代法計(jì)算時(shí)：令 fx=4x3-1.640x2+1.64x-0.410則 f'x=12x2-3.280x+1.64迭代公式為：xk+1=xk- fxf'x令x0=0.2進(jìn)行迭代，最后計(jì)算得氧氣分壓為：p=0.274901，則總壓為：P=3p+(1-2p)=1.274901二分法的計(jì)算算法：#include <iostream>#include<math.h>using namespace std;double a=0.2,b=0.3,e=0.001,n;double f(fl

4、oat c) n=4*c*c*c-1.6400000*c*c+1.64*c-0.410;return n;void main()double c,m,p,a=0.2,b=0.3;for(p=a;b-a>=0.001;)c=(a+b)/2;m=f(c);if(m<0)a=c;else if(m>0)b=c;elsep=c;break;p=(a+b)/2;cout<<" p="<<p<<endl;cout<<"n a="<<a<<endl; cout<<&q

5、uot;n b="<<b<<endl;牛頓迭代法的算法：#include<iostream>#include<string>#include<cmath>using namespace std;double f(double x)double m,n;m=4*x*x*x-1.640*x*x+1.64*x-0.410;n=12*x*x-3.280*x+1.64;return x-m/n;void newton(double x,double d)double a=x;double b=f(a);int k; /記錄循環(huán)的次數(shù)f

6、or(k=1;fabs(a-b)>d;k+)a=b;b=f(a);if(k>100)cout<<"迭代失敗，該函數(shù)可能不收斂!"<<endl;cout<<" a="<<a<<endl;cout<<"n b="<<b<<endl; cout<<"n k="<<k<<endl;return;int main()cout<<"請(qǐng)輸入初始值x0和要求得結(jié)果的精

7、度："double x,d;cin>>x>>d;newton(x,d);return 0;四 數(shù)值實(shí)驗(yàn)過(guò)程 二分法通過(guò)VC6.0程序運(yùn)行的結(jié)果如下圖所示： 牛頓迭代法通過(guò)VC6.0的計(jì)算結(jié)果如下：五 相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析由二分法程序輸出結(jié)果可得:a=0.274219,b=0.275。所以實(shí)驗(yàn)誤差：|x*-xk|（b-a）/2=0.0004<0.001 滿足實(shí)驗(yàn)要求。由牛頓迭代法程序輸出結(jié)果為a=0.274918，b=0.274901。實(shí)驗(yàn)誤差： |x*-xk|<|a-b|=0.000017<0.001， 滿足實(shí)驗(yàn)要求。二 拉格朗日插值法

8、一 問(wèn)題簡(jiǎn)介 在化學(xué)實(shí)驗(yàn)中，通常測(cè)得的是一批離散數(shù)據(jù)，需要從這批有限的測(cè)量數(shù)據(jù)中得出一個(gè)函數(shù)關(guān)系，進(jìn)而來(lái)求解任意的函數(shù)值，如：實(shí)驗(yàn)測(cè)得某物質(zhì)在20下，其粘度（Pa·s）與水溶液濃度c（重量%）的關(guān)系如下表。試用拉格朗日5次插值計(jì)算粘度在2.0×10-3和5.5×10-3時(shí)所對(duì)應(yīng)的濃度。要求精確到小數(shù)點(diǎn)后四位數(shù)。粘度（Pa·s）c（重量%）粘度（Pa·s）c（重量%）1.005×10-303.652×10-3401.776×10-3204.621×10-3452.480×10-3305.921&#

9、215;10-350二 數(shù)學(xué)模型在一定溫度下，物質(zhì)的水溶液粘度與其濃度是呈正相關(guān)的，表中給出6個(gè)點(diǎn)的函數(shù)值，所以可以進(jìn)行拉格朗日5次插值計(jì)算，其計(jì)算公式為：Lnx=k=0nykn+1(x)(x-xk)n+1'(xk)其中x即為粘度，Lnx即為所求值溶液濃度。三 算法選擇與算法過(guò)程計(jì)算算法：#include<iostream> using namespace std; const int N=1000; float Lagrange(float arrX,float arrY,int n,float x) float yResult=0.0; float LValueN; i

10、nt k,m; float temp1,temp2; for(k=0;k<n;k+) temp1=1.0; temp2=1.0; for(m=0;m<n;m+) if(m=k) continue; temp1 *= (x-arrXm); temp2 *= (arrXk-arrXm); LValuek=temp1/temp2; for(int i=0;i<n;i+) yResult += arrYi*LValuei; return yResult; int main() float arrXN,arrYN; int num; cout<<"輸入插值節(jié)點(diǎn)的個(gè)

11、數(shù)(小于"<<N<<"個(gè)): " cin>>num; cout<<"n-接下來(lái)輸入這些插值節(jié)點(diǎn)(先輸入X 再輸入對(duì)應(yīng)的Y)-n" for(int i=0;i<num;i+) cout<<"第"<<i+1<<"個(gè)節(jié)點(diǎn)的X值: " cin>>arrXi; cout<<"第"<<i+1<<"個(gè)節(jié)點(diǎn)的Y值: " cin>>ar

12、rYi; float X; cout<<"n-請(qǐng)輸入第一個(gè)待求解的插值節(jié)點(diǎn)的X值-n" cin>>X; float Res1 = Lagrange(arrX,arrY,num,X); cout<<"n-第一個(gè)點(diǎn)插值結(jié)果為: "<<Res1<<endl;cout<<"n-請(qǐng)輸入第二個(gè)待求解的插值節(jié)點(diǎn)的X值-n" cin>>X; float Res2 = Lagrange(arrX,arrY,num,X); cout<<"n-第二個(gè)點(diǎn)

13、插值結(jié)果為: "<<Res2<<endl; return 0; 四 數(shù)值實(shí)驗(yàn)過(guò)程五 相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析在實(shí)驗(yàn)過(guò)程中，拉格朗日高次差值對(duì)數(shù)據(jù)的誤差很敏感，由于舍入誤差，在計(jì)算中這種誤差很可能被放大，因此在實(shí)際工作中，節(jié)點(diǎn)數(shù)一般不超過(guò)七個(gè)。三 最小二乘法一 問(wèn)題簡(jiǎn)介用紫外吸收分光光度計(jì)測(cè)定牛血清蛋白的紫外吸收標(biāo)準(zhǔn)曲線時(shí)，先配制1.0mg/mL的牛血清蛋白溶液，然后分別稀釋為0.125mg/mL, 0.25mg/mL, 0.375mg/mL,0.5mg/mL，0.625mg/mL，0.75mg/mL，0.875mg/mL的溶液，分別測(cè)定不同濃度溶液在278n

14、m處的紫外吸收，實(shí)驗(yàn)數(shù)據(jù)記錄如下表，試?yán)米钚《朔ㄓ?jì)算出濃度c與吸光度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ì)及入射光的波長(zhǎng)有關(guān)，c為吸光物質(zhì)的濃度 ，b為吸收層厚度。當(dāng)保證入射光波長(zhǎng)與吸收層厚度不變時(shí)，吸光度A與濃度c是呈線性相關(guān)的。三 算法選擇與算法過(guò)程 選用最小二乘法做線性擬合。計(jì)算算法：#include<stdio.h>#i

15、nclude<malloc.h>double f(double a, int b)int i;double k = 1;for (i=0; i<b; i+)k *= a;return k;int main(void)int i, j, n, m;double a22, d2=0, e2, t; printf("請(qǐng)輸入Xi和Yi的個(gè)數(shù)：n");scanf("%d", &n);double *b = (double *) malloc(sizeof(double)*n);double *c = (double *) malloc(s

16、izeof(double)*n); printf("請(qǐng)輸入Xi：n");for (i=0; i<n; i+)scanf("%lf", &bi); printf("請(qǐng)輸入Yi：n");for (j=0; j<n; j+)scanf("%lf", &cj);for (i=0; i<2; i+)for (j=0; j<2; j+)aij = 0;for (m=0; m<n; m+)aij += f(bm,i*2+j);a11 = a10;a10 = a01;for (i=0;

17、 i<n; i+)d0 += ci; d1 += bi*ci;for(i=0; i<2; +i)t = a0i; a0i = a1i; a1i = t;t = d0; d0 = d1; d1 = t; a11 -= a01*a10/a00;d1 -= d0*a10/a00;a10 = 0;if(a11 = 0 && d1 !=0)printf("無(wú)解。");else e1 = d1/a11;e0 = (d0 - e1*a01)/a00; printf("方程的系數(shù)為：n"); for (j=0; j<2; j+)prin

18、tf("A%d = %10.6lfn",j+1, ej); printf("方程為：n"); printf("P(x) = %10.6lf + %10.6lfxn", e0, e1); free(b); free(c);return 0;四 數(shù)值實(shí)驗(yàn)過(guò)程五 相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析 由表中數(shù)據(jù)畫(huà)出的A-c的圖像如下圖所示： 在分析化學(xué)工作中，物質(zhì)的濃度與溶液的吸光度、電位或其它物理性質(zhì)存在線性關(guān)系，因此在測(cè)量中經(jīng)常要計(jì)算標(biāo)準(zhǔn)工作曲線法，而最小二乘法是在求標(biāo)準(zhǔn)曲線時(shí)的一種方便有效的方法。四 復(fù)合辛普森求積公式一 問(wèn)題簡(jiǎn)介 用脈沖法測(cè)定

19、某反應(yīng)器的停留時(shí)間分布，由實(shí)驗(yàn)數(shù)據(jù)計(jì)算壽命分布密度與時(shí)間t的關(guān)系如下表t10152025303540E(t)×10-301.929.020.9531.1032.2030.60求分布密度積分：1040E(t)dt二 數(shù)學(xué)模型 根據(jù)所給數(shù)據(jù)，小區(qū)間分為4等分，步長(zhǎng)h=10，n=4，用復(fù)合辛普森求積公式計(jì)算。Sn=h6fa+fb+4k=0n-1fxk+12+2k=1n-1fxkf(a)=f(10)=0 f(b)=f(40)=30.60×10-31040E(t)dt1060+30.6+41.92+20.95+32.2+29.0+31.1×10-3=0.5518三 算法過(guò)程

20、由于本題中的函數(shù)值都為散點(diǎn)，令m=2k，此時(shí)n=7，h=5則：Sn=h3fa+fb+4m=1,3,n-1fxm+2m=2,4,n-1fxm計(jì)算算法：#include<iostream>using namespace std;double Xinpusen(double a,double b,int n)int m,t,k;double 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=0;k<n/2;k+)m=2*k+1;s+=4*ym;for(k=1;k<n/2;

21、k+)t=2*k;s+=2*yt;return s*0.001*h/3;int main()double a=10;double b=40;int n=6; cout<<"由復(fù)化辛普森公式求的結(jié)果："<<Xinpusen(a,b,n)<<endl;return 0; 四 數(shù)值實(shí)驗(yàn)過(guò)程五 相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析 本實(shí)驗(yàn)通過(guò)10到40分鐘的反應(yīng)壽命分布密度數(shù)據(jù)來(lái)計(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

22、分鐘內(nèi)的反應(yīng)停留時(shí)間占總反應(yīng)時(shí)間的比率為55.18%。五 改進(jìn)歐拉公式一 問(wèn)題簡(jiǎn)介已知在管式反應(yīng)器中進(jìn)行液相反應(yīng)AR+S，反應(yīng)為吸熱反應(yīng)，反應(yīng)管外油浴溫度為340，假定已知管內(nèi)溫度與轉(zhuǎn)化率的關(guān)系為：dtdxA=-65.0-15.58(t-tc)/(k(1-xA)速度常數(shù)為：k=1.17×1017exp(-22400/T)tc=反應(yīng)器外壁溫度。如反應(yīng)器入口溫度為340，反應(yīng)器出口轉(zhuǎn)化率為90%，試用改進(jìn)歐拉公式求反應(yīng)器出口溫度。二 數(shù)學(xué)模型由題意知為常微分方程的初值問(wèn)題，假定tc=340(管外油浴溫度)，方程形式為：dtdxA=-65.0-15.58(t-340)/(k(1-xA)k=

23、1.17×1017exp(-22400/(t+273.2)xA=0，t0=340使用改進(jìn)歐拉法計(jì)算，設(shè)自變量xA的步長(zhǎng)h=0.05，經(jīng)過(guò)18步，可求出xA=0.9時(shí)的出口溫度。其數(shù)學(xué)模型為：yp=yn+h(yn-2xnyn)yc=yn+h(yp-2xn+1yp)yn+1=12(yp+yc)三 算法選擇與算法過(guò)程使用改進(jìn)歐拉法計(jì)算，設(shè)自變量xA的步長(zhǎng)h=0.05，經(jīng)過(guò)18步，可求出xA=0.9時(shí)的出口溫度。程序分為兩部分：(a)定義方程右邊函數(shù)；(b)用改進(jìn)歐拉法計(jì)算。計(jì)算算法：#include<iostream.h>#include<cmath>#define

24、 N 18void ModEuler(float (*f1)(float,float),float x0,float y0,float xn,int n)int i;float yp,yc,x=x0,y=y0,h=(xn-x0)/n;cout<<"x0="<<x<<'t'<<"y0"<<y<<endl;for(i=1;i<=n;i+)yp=y+h*f1(x,y);x=x0+i*h;yc=y+h*f1(x,yp);y=(float)(yp+yc)/2.0;cout

25、<<"x"<<i<<"="<<x<<" y"<<i<<"="<<y<<endl;void main()float xn=0.9,x0=0.0,y0=340;float f1(float ,float);ModEuler(f1,x0,y0,xn,N);float f1(float x,float y)float m=-22400/(y+273.2);float k=1.17*pow(10,17)*exp(m);

26、return -65.50-15.58*(y-340)/(k*(1-x);四 數(shù)值實(shí)驗(yàn)過(guò)程五 相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析 化學(xué)中最常見(jiàn)的常微分方程問(wèn)題為涉及傳熱和擴(kuò)散的問(wèn)題，如間歇反應(yīng)器的計(jì)算、活塞流反應(yīng)器的計(jì)算、定態(tài)一維熱傳導(dǎo)問(wèn)題、固定床反應(yīng)器的分散模型等問(wèn)題，都可以應(yīng)用本實(shí)驗(yàn)中的改進(jìn)歐拉公式進(jìn)行求解。六 列主元消去法一 問(wèn)題簡(jiǎn)介 用酸堿滴定法分析混合堿中Na2CO3、NaHCO3和NaOH的含量時(shí)，稱取混合堿0.2000g溶于適量水中，用0.1mol/L的標(biāo)準(zhǔn)HCl滴定，用酚酞指示劑由紅色轉(zhuǎn)無(wú)色，消耗HCl8.01mol，再用甲基橙指示繼續(xù)滴定到終點(diǎn)，共消耗HCl30.0mL，試計(jì)算混合堿

27、中Na2CO3、NaHCO3和NaOH的含量。二 數(shù)學(xué)模型設(shè)混合堿中Na2CO3、NaHCO3和NaOH的質(zhì)量（mg）分別為x1、x2、x3，查得Na2CO3、NaHCO3和NaOH的分子量分別為100、80、40。由題意可列出如下線性方程組：x1+x2+x3=200x1100+x340=0.1×82x1100+x280+x340=0.1×30可令A(yù)=1111/10001/402/1001/401/40 x=x1x2x3 b=2000.83則： Ax=b三 算法選擇與算法過(guò)程選用列主元消去法解上述線性方程組。A=1110.0100.0250.020.01250.025計(jì)算算

28、法：#include<iostream>#include<cstdio>#include<iomanip>using namespace std;#define maxn 50int n;double amaxnmaxn;double bmaxn;/double xmaxn;/void read() cout<<"請(qǐng)輸入系數(shù)矩陣規(guī)模n:= " cin>>n; cout<<"|請(qǐng)輸入系數(shù)矩陣：n" for(int i=1;i<=n;i+) for(int j=1;j<=n;

29、j+) cin>>aij; cout<<"|請(qǐng)輸入b矩陣：n" for(i=1;i<=n;i+) cin>>bi; void Print() cout<<"|-n" cout<<"|結(jié)果為：n" for(int i=1;i<=n;i+) cout<<"x"<<i<<"="<<xi<<endl; cout<<"|-nn"void Lie

30、ZhuXiaoYuan() for(int k=1;k<n;k+) double ab_max=-1; int max_ik; for(int i=k;i<=n;i+) if(abs(aik)>ab_max) ab_max=abs(aik); max_ik=i; if(max_ik!=k) double temp; for(int j=1;j<=n;j+) temp=amax_ikj; amax_ikj=akj; akj=temp; temp=bmax_ik; bmax_ik=bk; bk=temp; for(i=k+1;i<=n;i+) aik/=akk; fo

31、r(int j=k+1;j<=n;j+) aij-=aik*akj; bi-=aik*bk; if(k<n-1)continue; elsexn=bn/ann;for(int i=n-1;i>0;i-)xi=bi;for(int j=i+1;j<=n;j+)xi-=aij*xj;xi/=aii; Print();int main() while(1) read(); LieZhuXiaoYuan(); 四 數(shù)值實(shí)驗(yàn)過(guò)程五 相關(guān)數(shù)值分析和實(shí)際應(yīng)用分析 高斯消去法既容易以理解又很簡(jiǎn)單，但經(jīng)常得不到正確結(jié)果，原因是消元過(guò)程是按方程組的自然順序進(jìn)行的。如果消去元數(shù)值太小，會(huì)產(chǎn)生較大的舍入誤差，得出完全錯(cuò)誤的結(jié)果。如果消去元是零會(huì)導(dǎo)致消元過(guò)程無(wú)法進(jìn)行。為了解決這些問(wèn)題，就要選用列主元消去法。在根據(jù)光度確定多組分溶液中各組分濃度時(shí)也常常用到列主元消去法。七 不動(dòng)點(diǎn)迭代法一 問(wèn)題簡(jiǎn)介 在9.33atm，302K時(shí)，容器中充2mol氨氣，試求該容器的容積，要求精度為10-2。（氨氣的范德華常數(shù)：a=4.17atm·L