数值分析报告上机C版

《数值分析报告上机C版》由会员分享，可在线阅读，更多相关《数值分析报告上机C版（30页珍藏版）》请在装配图网上搜索。

1、word数值分析论文NUMERICAL ANALYSIS(THESIS)题目数值分析论文学生指导教师鸿雁学院专业班级2016年12月目录一 二分法与牛顿迭代法1二 拉格朗日插值法6三 最小二乘法10四 复合辛普森求积公式15五 改良欧拉公式18六 列主元消去法21七 不动点迭代法2628 / 30一 二分法与牛顿迭代法一问题简介在有气体参加的恒容反响体系里，反响体系的总压会随着反响的进展而改变。通过反响的平衡常数可以求出反响后气体的总量，而气体分子的分压与其所占比率是成正比的。这类问题中关键是要求出反响后各分子的分压，各分子分压之和即为总压，其实质即为求方程的根。例如：在298K下，化学反响2

2、OF2=O2+2F2的平衡常数为0.410atm，如在298K下将OF2通入容器，当t=0时为1atm，问最后总压是多少？计算精度为10-3。二 数学模型假设是理想气体，可知 2OF2=O2+2F2设氧的分压为p，平衡时有 1-2p p p平衡时有 整理得 4p3-1640p2+164p-0410=0函数关系式为 三 算法选择与算法过程由计算得 f(0.2)=因此，有根区间为0.2，0.3用求单根的二分法计算时：可令a=0.2，b=0.3；最后计算609atm，总压为：609atm用牛顿迭代法计算时：令如此迭代公式为：令进展迭代，最后计算得氧气分压为：p=0.274901，如此总压为：P=二分

3、法的计算算法：#include #includeusing namespace std;double a=0.2,b=0.3,e=0.001,n;double f(float 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(m0)b=c;elsep=c;break;p=(a+b)/2;cout p=pendl;coutn a=aendl; coutn b=bendl;牛顿迭代法的算法：#inclu

4、de#include#includeusing 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; /记录循环的次数for(k=1;fabs(a-b)d;k+)a=b;b=f(a);if(k100)cout迭代失败，该函数可能不收敛!endl;cout a=aendl;coutn b=bendl; cou

5、tn k=kendl;return;int main()coutxd;newton(x,d);return 0;四 数值实验过程二分法通过VC6.0程序运行的结果如如下图所示：牛顿迭代法通过VC6.0的计算结果如下：五 相关数值分析和实际应用分析由二分法程序输出结果可得:a=0.274219,b=0.275。所以实验误差：|x*-xk|b-a/2 满足实验要求。由牛顿迭代法程序输出结果为a=0.274918，b=0.274901。实验误差： |x*-xk|a-b|=0.0000170.001， 满足实验要求。二 拉格朗日插值法一 问题简介 在化学实验中，通常测得的是一批离散数据，需要从这批有限

6、的测量数据中得出一个函数关系，进而来求解任意的函数值，如：实验测得某物质在20下，其粘度Pas与水溶液浓度c重量%10-310-3时所对应的浓度。要求准确到小数点后四位数。粘度Pasc重量%粘度Pasc重量%10-3010-34010-32010-34510-33010-350二 数学模型在一定温度下，物质的水溶液粘度与其浓度是呈正相关的，表中给出6个点的函数值，所以可以进展拉格朗日5次插值计算，其计算公式为：其中x即为粘度，即为所求值溶液浓度。三 算法选择与算法过程计算算法：#include using namespace std; const int N=1000; float Lagra

7、nge(float arrX,float arrY,int n,float x) float yResult=0.0; float LValueN; int k,m; float temp1,temp2; for(k=0;kn;k+) temp1=1.0; temp2=1.0; for(m=0;mn;m+) if(m=k) continue; temp1 *= (x-arrXm); temp2 *= (arrXk-arrXm); LValuek=temp1/temp2; for(int i=0;in;i+) yResult += arrYi*LValuei; return yResult; i

8、nt main() float arrXN,arrYN; int num; cout输入插值节点的个数(小于Nnum; coutn-接下来输入这些插值节点(先输入X 再输入对应的Y)-n; for(int i=0;inum;i+) cout第i+1arrXi; cout第i+1arrYi; float X; coutX; float Res1 = Lagrange(arrX,arrY,num,X); coutn-第一个点插值结果为: Res1endl;coutX; float Res2 = Lagrange(arrX,arrY,num,X); coutn-第二个点插值结果为: Res2endl

9、; return 0; 四 数值实验过程五 相关数值分析和实际应用分析在实验过程中，拉格朗日高次差值对数据的误差很敏感，由于舍入误差，在计算中这种误差很可能被放大，因此在实际工作中，节点数一般不超过七个。三 最小二乘法一 问题简介用紫外吸收分光光度计测定牛血清蛋白的紫外吸收标准曲线时，先配制mg/mL, mg/mL, mg/mL,的溶液，分别测定不同浓度溶液在278nm处的紫外吸收，实验数据记录如下表，试利用最小二乘法计算出浓度c与吸光度A之间的关系。cAcA1二 数学模型物质的紫外吸收吸光度与它的浓度满足朗伯比尔定律,即：式中A为吸光度，K为摩尔吸收系数，它与吸收物质的性质与入射光的波长有关

10、，c为吸光物质的浓度 ，b为吸收层厚度。当保证入射光波长与吸收层厚度不变时，吸光度A与浓度c是呈线性相关的。三 算法选择与算法过程 选用最小二乘法做线性拟合。计算算法：#include#includedouble f(double a, int b)int i;double k = 1;for (i=0; ib; i+)k *= a;return k;int main(void)int i, j, n, m;double a22, d2=0, e2, t; printf(请输入Xi和Yi的个数：n);scanf(%d, &n);double *b = (double *) malloc(siz

11、eof(double)*n);double *c = (double *) malloc(sizeof(double)*n); printf(请输入Xi：n);for (i=0; in; i+)scanf(%lf, &bi); printf(请输入Yi：n);for (j=0; jn; j+)scanf(%lf, &cj);for (i=0; i2; i+)for (j=0; j2; j+)aij = 0;for (m=0; mn; m+)aij += f(bm,i*2+j);a11 = a10;a10 = a01;for (i=0; in; i+)d0 += ci; d1 += bi*ci;

12、for(i=0; i2; +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(无解。);else e1 = d1/a11;e0 = (d0 - e1*a01)/a00; printf(方程的系数为：n); for (j=0; j2; j+)printf(A%d = %10.6lfn,j+1, ej); printf(方程为：n); printf(P(x) = %10.6lf + %10.6lfx

13、n, e0, e1); free(b); free(c);return 0;四 数值实验过程五 相关数值分析和实际应用分析 由表中数据画出的A-c的图像如如下图所示：在分析化学工作中，物质的浓度与溶液的吸光度、电位或其它物理性质存在线性关系，因此在测量中经常要计算标准工作曲线法，而最小二乘法是在求标准曲线时的一种方便有效的方法。四 复合辛普森求积公式一 问题简介 用脉冲法测定某反响器的停留时间分布，由实验数据计算寿命分布密度与时间t的关系如下表t10152025303540E(t)10-30求分布密度积分：二 数学模型 根据所给数据，小区间分为4等分，步长h=10，n=4，用复合辛普森求积公式

14、计算。f(a)=f(10)=0 f(b)=f(40)10-3三 算法过程由于此题中的函数值都为散点，令m=2k，此时n=7，h=5如此：计算算法：#includeusing 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;kn/2;k+)m=2*k+1;s+=4*ym;for(k=1;kn/2;k+)t=2*k;s+=2*yt;return s

15、*0.001*h/3;int main()double a=10;double b=40;int n=6; cout由复化辛普森公式求的结果：Xinpusen(a,b,n)endl;return 0; 四 数值实验过程五 相关数值分析和实际应用分析 本实验通过10到40分钟的反响寿命分布密度数据来计算反响器的停留时间分布，利用复合辛普森求积法对10到40分钟的寿命分布密度进展积分，即得10到40分钟反响停留时间占总反响时间的比例。实验计算结果为0.5518，即10到40分钟的反响停留时间占总反响时间的比率为55.18%。五 改良欧拉公式一 问题简介在管式反响器中进展液相反响AR+S，反响为吸热

16、反响，反响管外油浴温度为340，假定管温度与转化率的关系为：1017exp(-22400/T)tc=反响器外壁温度。如反响器入口温度为340，反响器出口转化率为90%，试用改良欧拉公式求反响器出口温度。二 数学模型由题意知为常微分方程的初值问题，假定tc=340(管外油浴温度)，方程形式为：1017exp(-22400/(t+273.2)xA=0，t0=340使用改良欧拉法计算，设自变量xA的步长h=0.05，经过18步，可求出xA=0.9时的出口温度。其数学模型为：三 算法选择与算法过程使用改良欧拉法计算，设自变量xA的步长h=0.05，经过18步，可求出xA=0.9时的出口温度。程序分为两

17、局部：(a)定义方程右边函数；(b)用改良欧拉法计算。计算算法：#include#include#define 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;coutx0=xty0yendl;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;coutxi=x yi=yendl;void main()float

18、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);return -65.50-15.58*(y-340)/(k*(1-x);四 数值实验过程五 相关数值分析和实际应用分析 化学中最常见的常微分方程问题为涉与传热和扩散的问题，如间歇反响器的计算、活塞流反响器的计算、定态一维热传导问题、固定床反响器的分散模型等问题，都可以应用本实验中的改良欧拉公式进展求

19、解。六 列主元消去法一 问题简介 用酸碱滴定法分析混合碱中Na2CO3、NaHCO3和NaOH的含量时，称取混合碱0.2000g溶于适量水中，用0.1mol/L的标准HCl滴定，用酚酞指示剂由红色转无色，消耗HCl8.01mol，再用甲基橙指示继续滴定到终点，共消耗HCl30.0mL，试计算混合碱中Na2CO3、NaHCO3和NaOH的含量。二 数学模型设混合碱中Na2CO3、NaHCO3和NaOH的质量mg分别为x1、x2、x3，查得Na2CO3、NaHCO3和NaOH的分子量分别为100、80、40。由题意可列出如下线性方程组：可令如此： Ax=b三 算法选择与算法过程选用列主元消去法解上

20、述线性方程组。计算算法：#include#include#includeusing namespace std;#define maxn 50int n;double amaxnmaxn;double bmaxn;/double xmaxn;/void read() coutn; cout|请输入系数矩阵：n; for(int i=1;i=n;i+) for(int j=1;jaij; cout|请输入b矩阵：n; for(i=1;ibi; void Print() cout|-n; cout|结果为：n; for(int i=1;i=n;i+) coutxi=xiendl; cout|-nn

21、;void LieZhuXiaoYuan() for(int k=1;kn;k+) double ab_max=-1; int max_ik; for(int i=k;iab_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; for(int j=k+1;j=n;j+) aij-=aik*

22、akj; bi-=aik*bk; if(k0;i-)xi=bi;for(int j=i+1;j=n;j+)xi-=aij*xj;xi/=aii; Print();int main() while(1) read(); LieZhuXiaoYuan(); 四 数值实验过程五 相关数值分析和实际应用分析高斯消去法既容易以理解又很简单，但经常得不到正确结果，原因是消元过程是按方程组的自然顺序进展的。如果消去元数值太小，会产生较大的舍入误差，得出完全错误的结果。如果消去元是零会导致消元过程无法进展。为了解决这些问题，就要选用列主元消去法。在根据光度确定多组分溶液中各组分浓度时也常常用到列主元消去法。七

23、 不动点迭代法一 问题简介 在9.33atm，302K时，容器中充2mol氨气，试求该容器的容积，要求精度为10-2L2/mol2atmL2/mol2二 数学模型德华方程：改写为迭代公式并代入数据：初值V0用理想气体状态方程估计：用初值V0=5.2833L带入迭代公式求解。构造迭代方程： 初值： 三 算法选择与算法过程选用不动点迭代法进展计算计算算法：#include#include#includeusing namespace std;float f(float x)float a=9.33+16.68/(x*x);return 49.2948/a+0.0742;void Diedai(float x,float d)float a=x;float b=f(a);int k; /记录循环的次数for(k=1;a-bd;k+)a=b;b=f(a);coutxk=a100)cout迭代失败！可能是函数不收敛endl;return ;coutb=bendl;cout迭代次数为kendl;return;int main()coutxd;Diedai(x,d);return 0;四 数值实验过程实验计算结果如下列图：由实验结果得：V=5.0064，迭代次数为3五 相关数值分析和实际应用分析收敛性分析： 如此迭代法是收敛的。