数值分析上机(C++版)

数值分析论文NUMERICAL ANALYSIS(THESIS)题 目数值分析论文学生姓名指导教师张鸿雁学 院专业班级2016年12月目录一 二分法与牛顿迭代法1二 拉格朗日插值法 6三 最小二乘法10四 复合辛普森求积公式15五 改进欧拉公式18六 列主元消去法21七 不动点迭代法26一 二分法与牛顿迭代法一 问题简介在有气体参加的恒容反应体系里，反应体系的总压会随着反应的进行而改变。通过反应的平衡常数可以求出反应后气体的总量，而气体分子的分压与其所占比率是成正比的。这类问题中关键是要求出反应后各分子的分压，各分子分压之和即为总压，其实质即为求方程的根。例如：在298K下，化学反应2OF2=O2+2F2的平衡常数为0.410atm，如在298K下将OF2通入容器，当t=0时为1atm，问最后总压是多少？计算精度为10-3。二 数学模型假设是理想气体，可知 2OF2=O2+2F2设氧的分压为p，平衡时有 1-2p p p平衡时有 4p3(1-2p)2=0.410整理得 4p3-1.640p2+1.64p-0.410=0函数关系式为 fp=4p3-1.640p2+1.64p-0.410=0三 算法选择与算法过程由计算得 f(0.2)=0.1156，f(0.3)=0.0424因此，有根区间为0.2，0.3用求单根的二分法计算时：可令a=0.2，b=0.3；最后计算得氧气分压：p=0.274609atm，总压为：P=3p+(1-2p)=1.274609atm用牛顿迭代法计算时：令 fx=4x3-1.640x2+1.64x-0.410则 fx=12x2-3.280x+1.64迭代公式为：xk+1=xk- fxfx令x0=0.2进行迭代，最后计算得氧气分压为：p=0.274901，则总压为：P=3p+(1-2p)=1.274901二分法的计算算法：#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;牛顿迭代法的算法：#include#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; coutn 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=0.00040.001 满足实验要求。由牛顿迭代法程序输出结果为a=0.274918，b=0.274901。实验误差： |x*-xk|a-b|=0.0000170.001， 满足实验要求。二 拉格朗日插值法一 问题简介 在化学实验中，通常测得的是一批离散数据，需要从这批有限的测量数据中得出一个函数关系，进而来求解任意的函数值，如：实验测得某物质在20下，其粘度（Pas）与水溶液浓度c（重量%）的关系如下表。试用拉格朗日5次插值计算粘度在2.010-3和5.510-3时所对应的浓度。要求精确到小数点后四位数。粘度（Pas）c（重量%）粘度（Pas）c（重量%）1.00510-303.65210-3401.77610-3204.62110-3452.48010-3305.92110-350二 数学模型在一定温度下，物质的水溶液粘度与其浓度是呈正相关的，表中给出6个点的函数值，所以可以进行拉格朗日5次插值计算，其计算公式为：Lnx=k=0nykn+1(x)(x-xk)n+1(xk)其中x即为粘度，Lnx即为所求值溶液浓度。三 算法选择与算法过程计算算法：#include using namespace std; const int N=1000; float Lagrange(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; int 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; return 0; 四 数值实验过程五 相关数值分析和实际应用分析在实验过程中，拉格朗日高次差值对数据的误差很敏感，由于舍入误差，在计算中这种误差很可能被放大，因此在实际工作中，节点数一般不超过七个。三 最小二乘法一 问题简介用紫外吸收分光光度计测定牛血清蛋白的紫外吸收标准曲线时，先配制1.0mg/mL的牛血清蛋白溶液，然后分别稀释为0.125mg/mL, 0.25mg/mL, 0.375mg/mL,0.5mg/mL，0.625mg/mL，0.75mg/mL，0.875mg/mL的溶液，分别测定不同浓度溶液在278nm处的紫外吸收，实验数据记录如下表，试利用最小二乘法计算出浓度c与吸光度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是呈线性相关的。三 算法选择与算法过程 选用最小二乘法做线性拟合。计算算法：#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(sizeof(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;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.6lfxn, e0, e1); free(b); free(c);return 0;四 数值实验过程五 相关数值分析和实际应用分析 由表中数据画出的A-c的图像如下图所示： 在分析化学工作中，物质的浓度与溶液的吸光度、电位或其它物理性质存在线性关系，因此在测量中经常要计算标准工作曲线法，而最小二乘法是在求标准曲线时的一种方便有效的方法。四 复合辛普森求积公式一 问题简介 用脉冲法测定某反应器的停留时间分布，由实验数据计算寿命分布密度与时间t的关系如下表t10152025303540E(t)10-301.929.020.9531.1032.2030.60求分布密度积分：1040E(t)dt二 数学模型 根据所给数据，小区间分为4等分，步长h=10，n=4，用复合辛普森求积公式计算。Sn=h6fa+fb+4k=0n-1fxk+12+2k=1n-1fxkf(a)=f(10)=0 f(b)=f(40)=30.6010-31040E(t)dt1060+30.6+41.92+20.95+32.2+29.0+31.110-3=0.5518三 算法过程由于本题中的函数值都为散点，令m=2k，此时n=7，h=5则：Sn=h3fa+fb+4m=1,3,n-1fxm+2m=2,4,n-1fxm计算算法：#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*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，反应为吸热反应，反应管外油浴温度为340，假定已知管内温度与转化率的关系为：dtdxA=-65.0-15.58(t-tc)/(k(1-xA)速度常数为：k=1.171017exp(-22400/T)tc=反应器外壁温度。如反应器入口温度为340，反应器出口转化率为90%，试用改进欧拉公式求反应器出口温度。二 数学模型由题意知为常微分方程的初值问题，假定tc=340(管外油浴温度)，方程形式为：dtdxA=-65.0-15.58(t-340)/(k(1-xA)k=1.171017exp(-22400/(t+273.2)xA=0，t0=340使用改进欧拉法计算，设自变量xA的步长h=0.05，经过18步，可求出xA=0.9时的出口温度。其数学模型为：yp=yn+h(yn-2xnyn)yc=yn+h(yp-2xn+1yp)yn+1=12(yp+yc)三 算法选择与算法过程使用改进欧拉法计算，设自变量xA的步长h=0.05，经过18步，可求出xA=0.9时的出口温度。程序分为两部分：(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 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);四 数值实验过程五 相关数值分析和实际应用分析 化学中最常见的常微分方程问题为涉及传热和扩散的问题，如间歇反应器的计算、活塞流反应器的计算、定态一维热传导问题、固定床反应器的分散模型等问题，都可以应用本实验中的改进欧拉公式进行求解。六 列主元消去法一 问题简介 用酸碱滴定法分析混合碱中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。由题意可列出如下线性方程组：x1+x2+x3=200x1100+x340=0.182x1100+x280+x340=0.130可令A=1111/10001/402/1001/401/40 x=x1x2x3 b=2000.83则： Ax=b三 算法选择与算法过程选用列主元消去法解上述线性方程组。A=1110.0100.0250.020.01250.025计算算法：#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;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*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(); 四 数值实验过程五 相关数值分析和实际应用分析 高斯消去法既容易以理解又很简单，但经常得不到正确结果，原因是消元过程是按方程组的自然顺序进行的。如果消去元数值太小，会产生较大的舍入误差，得出完全错误的结果。如果消去元是零会导致消元过程无法进行。为了解决这些问题，就要选用列主元消去法。在根据光