浙江大学计算方法大作业第十二题

题目：某疾病发生率y‰和年龄段x（每五年为一段，例如0～5岁为第一段，6～10岁为第二段，……）之间有形如=bxyae的关系。试根据观察所得的如下数据表，用最小二乘法确定式中的参数a和b，并计算相应的均方误差和最大偏差。x12345678910y0.8982.383.071.842.021.942.222.774.024.76x111213141516171819y5.466.5310.916.522.535.750.661.681.8思路步骤：解线性方程组法1.根据x和y的直接关系bxyae=可得lnlnyabx=+；将求解的问题化归为求解lny与x之间的线性关系；2.令多元函数20(ln)miiiabxy=+-å对a，b的偏导等于0，求解线性方程组即可得到所求的a，b；3.将得到的a，b带回求解出均方误差与最大偏差。解非线性方程组1.令多元函数20()imbxiiaey=-å对a，b的偏导等于0，得到关于a，b的两个非线性方程组，(,)0,(,)0fabgab==，利用牛顿迭代法求解；2.取一组00,ab；3.令方程组：0000000000000000(,)()()0,(,)()()0;aaaabbbbaaaabbbbfffabaabbabgggabaabbab========抖+-+-=抖抖+-+-=抖4.求解上个步骤中得出的a，b；5.若00.00001aa-或者00.00001bb-则重复上述的两个步骤；6.当循环停止时得到所求的a，b。计算结果：a=0.681361b=0.230620均方误差为38.325470最大偏差为27.304749a=0.236882b=0.308978均方误差为8.655291最大偏差为5.339718结果分析：解决本题采用第一种两边取对数，然后求解线性方程组，但是这种方法求出的均方误差和最大偏差并不是在原指数函数的情况下满足的最优结果，所以均方误差和最大偏差的结果都不够理想。改进后的方法采用了牛顿迭代法，通过运用迭代法求解非线性方程组，求得的结果有比较好的收敛效果，同时使均方误差和最大偏差较小（均方误差和最大偏差分别为第一种方法的22.58%和19.56%），取得了比较理想的效果。源程序：#includestdio.h#includemath.hintmain(void){FILE*fp1,*fp2,*fp3;inti,j;doublex[19],y[19],M[2][2],N[2];doublea,b,lna;doubleA[2][2],B[2];doublea0,b0;doublesigma,max;if((fp1=fopen(数据输入.txt,r))==NULL){printf(Fileopenerror!\n);exit(0);}for(i=0;i19;i++)fscanf(fp1,%lf%lf,&x[i],&y[i]);if((fp2=fopen(数据输出1.txt,w))==NULL){printf(Fileopenerror!\n);exit(0);}if((fp3=fopen(数据输出2.txt,w))==NULL){printf(Fileopenerror!\n);exit(0);}sigma=max=0;M[0][0]=M[0][1]=M[1][0]=M[1][1]=N[0]=N[1]=0;for(i=0;i19;i++){M[0][0]=19;M[0][1]=M[0][1]+x[i];M[1][0]=M[1][0]+x[i];M[1][1]=M[1][1]+x[i]*x[i];N[0]=N[0]+1.0*log(y[i]);N[1]=N[1]+1.0*x[i]*log(y[i]);}lna=(N[0]*M[1][1]-N[1]*M[0][1])/(M[0][0]*M[1][1]-M[0][1]*M[1][0]);b=(N[1]*M[0][0]-N[0]*M[1][0])/(M[0][0]*M[1][1]-M[0][1]*M[1][0]);a=exp(lna);max=fabs(a*exp(b*x[1])-y[1]);for(i=0;i19;i++){sigma=sigma+pow(a*exp(b*x[i])-y[i],2);if(fabs(a*exp(b*x[i])-y[i])max)max=fabs(a*exp(b*x[i])-y[i]);}sigma=pow(sigma,0.5);fprintf(fp2,a=%lf\n,a);fprintf(fp2,b=%lf\n,b);fprintf(fp2,均方误差为%lf\n,sigma);fprintf(fp2,最大偏差为%lf\n,max);sigma=max=0;A[0][0]=A[0][1]=A[1][0]=A[1][1]=B[0]=B[1]=0;a0=0.24;b0=0.31;for(i=0;i19;i++){A[0][0]=A[0][0]+exp(2*b0*x[i]);A[0][1]=A[0][1]+2*a0*x[i]*exp(2*b0*x[i])-x[i]*y[i]*exp(b0*x[i]);B[0]=B[0]+(2*a0*x[i]*exp(2*b0*x[i])-x[i]*y[i]*exp(b0*x[i]))*b0+exp(2*b0*x[i])*a0+y[i]*exp(b0*x[i])-a0*exp(2*b0*x[i]);A[1][0]=A[1][0]+2*a0*x[i]*exp(2*b0*x[i])-x[i]*y[i]*exp(b0*x[i]);A[1][1]=A[1][1]+2*a0*a0*x[i]*x[i]*exp(2*b0*x[i])-a0*x[i]*x[i]*y[i]*exp(b0*x[i]);B[1]=B[1]+(2*a0*x[i]*exp(2*b0*x[i])-x[i]*y[i]*exp(b0*x[i]))*a0+(2*a0*a0*x[i]*x[i]*exp(2*b0*x[i])-a0*x[i]*x[i]*y[i]*exp(b0*x[i]))*b0+a0*x[i]*y[i]*exp(b0*x[i])-a0*a0*x[i]*exp(2*b0*x[i]);}a=(B[0]*A[1][1]-B[1]*A[0][1])/(A[0][0]*A[1][1]-A[0][1]*A[1][0]);b=(B[1]*A[0][0]-B[0]*A[1][0])/(A[0][0]*A[1][1]-A[0][1]*A[1][0]);while(fabs(a-a0)0.000000001||fabs(b-b0)0.000000001){a0=a;b0=b;A[0][0]=A[0][1]=A[1][0]=A[1][1]=B[0]=B[1]=0;for(i=0;i19;i++){A[0][0]=A[0][0]+exp(2*b0*x[i]);A[0][1]=A[0][1]+2*a0*x[i]*exp(2*b0*x[i])-x[i]*y[i]*exp(b0*x[i]);B[0]=B[0]+(2*a0*x[i]*exp(2*b0*x[i])-x[i]*y[i]*exp(b0*x[i]))*b0+exp(2*b0*x[i])*a0+y[i]*exp(b0*x[i])-a0*exp(2*b0*x[i]);A[1][0]=A[1][0]+2*a0*x[i]*exp(2*b0*x[i])-x[i]*y[i]*exp(b0*x[i]);A[1][1]=A[1][1]+2*a0*a0*x[i]*x[i]*exp(2*b0*x[i])-a0*x[i]*x[i]*y[i]*exp(b0*x[i]);B[1]=B[1]+(2*a0*x[i]*exp(2*b0*x[i])-x[i]*y[i]*exp(b0*x[i]))*a0+(2*a0*a0*x[i]*x[i]*exp(2*b0*x[i])-a0*x[i]*x[i]*y[i]*exp(b0*x[i]))*b0+a0*x[i]*y[i]*exp(b0*x[i])-a0*a0*x[i]*exp(2*b0*x[i]);}a=(B[0]*A[1][1]-B[1]*A[0][1])/(A[0][0]*A[1][1]-A[0][1]*A[1][0]);b=(B[1]*A[0][0]-B[0]*A[1][0])/(A[0][0]*A[1][1]-A[0][1]*A[1][0]);}max=fabs(a*exp(b*x[1])-y[1]);for(i=0;i19;i++){sigma=sigma+pow(a*exp(b*x[i])-y[i],2);if(fabs(a*exp(b*x[i])-y[i])max)max=fabs(a*exp(b*x[i])-y[i]);}sigma=pow(sigma,0.5);fprintf(fp3,a=%lf\n,a);fprintf(fp3,b=%lf\n,b);fprintf(fp3,均方误差为%lf\n,sigma);fprintf(fp3,最大偏差为%lf\n,max);}