科学计算 实习题三 曲线拟合的最小二乘法

1、 首先参照教材上的方法对方程 y=aebx 进行一定的变换，两边取对数可以得到
ln(y)=ln(a)+bx
令 y1=ln(y)，a1=ln(a)，得 y1=a1+bx 这一过程程序代码为
for(i=0;i<19;i++){ scanf("%lf",&y0[i]);
} for(j=0;j<19;j++){
e=(aeb-y1)2+(ae2b-y2)2+( ae3b-y3)2+…+( ae19b-y19)2
∂������
∂������
求出相应的方程组 f(a,b)=∂������ ，g(a,b)= ∂������，分别记为 f 与 g，再计算出法方
∂������ ∂������ ∂������ ∂������ 程组的各项偏导数∂������、∂������、∂������、∂������。分别记为 f1，f2，g1，g2。
程序代码为 a1=(s*p-t*m)/(p*p-m*q); a0=(s-p*a1)/m; b=a1; a=pow(e,a0);
4、 最后，按照题目要求求得这一组参数公式所对应的均方误差与最大偏差，这一过 程程序代码为 f=z=0; for(l=0;l<19;l++){ f=f+(y0[l]-a*pow(e,b*(x[l])))*(y0[l]-a*pow(e,b*(x[l]))); if(fabs((y0[l]-a*pow(e,b*(x[l]))))>z)z=fabs((y0[l]-a*pow(e,b*(x[l])))); } f0=sqrt(f);
printf("a=%lf,b=%lf\n 均方误差为%lf,最大偏差为%lf\n", a, b, f0, z); } 方法二、
#include<stdio.h> #include<math.h> int main(void) {
double a,b,e,m,f,f1,f2,g,g1,g2,c,d; double x[19],y[19]; int i,j;
f=f+2*(a*exp((i+1)*b)-y[i])*exp((i+1)*b); f1=f1+2*exp(2*(i+1)*b); f2=f2+4*(i+1)*a*exp(2*(i+1)*b)-2*(i+1)*y[i]*exp((i+1)*b); g=g+2*(a*exp((i+1)*b)-y[i])*(i+1)*a*exp((i+1)*b); g1=g1+4*(i+1)*a*exp(2*(i+1)*b)-2*(i+1)*y[i]*exp((i+1)*b); g2=g2+4*pow((i+1)*a,2)*exp(2*(i+1)*b)-2*pow((i+1),2)*a*y[i]*exp((i+1)*b); } c=f1*a+f2*b-f; d=g1*a+g2*b-g; b=(d-c*(g1/f1))/(g2-f2*(g1/f1)); a=(c-f2*b)/f1; 其中，j 为所需要的迭代次数。 6、 最后，根据得到的 a，b 的满意解，求出均方误差与最大偏差，这一过程程序代码 可参照方法一中的代码源程序。 三、 计算结果 1、
a1=(s*p-t*m)/(p*p-m*q); a0=(s-p*a1)/m; b=a1; a=pow(e,a0);
f=z=0; for(l=0;l<19;l++){
f=f+(y0[l]-a*pow(e,b*(x[l])))*(y0[l]-a*pow(e,b*(x[l]))); if(fabs((y0[l]-a*pow(e,b*(x[l]))))>z)z=fabs((y0[l]-a*pow(e,b*(x[l])))); } f0=sqrt(f);
科学计算上机实习题目三
------曲线拟合的最小二乘法
3120104128 杨超
一、 题目 某疾病发生率 y‰和年龄段 x（每五年为一段，例如 0~5 岁为第一段，6~10 岁为第二
段……）之间有形如 y=aebx 的关系。是根据观测得到的如下数据表，用最小二乘法
确定式中的参数 a 和 b，并计算相应的均方误差与最大偏差（注：本题应在方法应用 上需要经过一定努力，才能得到较理想结果）。
计算出 a，b 的值 （3）、以新得到的 a、b 的值为 a0，b0 带入（2）中公式进行迭代计算，直到得到 满意的解。计算过程的源程序为
scanf("%lf",&a); scanf("%lf",&b); f=f1=f2=g=g1=g2=0;
for(j=0;j<10;j++){ for(i=0;i<19;i++){
} c=f1*a+f2*b-f; d=g1*a+g2*b-g; b=(d-c*(g1/f1))/(g2-f2*(g1/f1)); a=(c-f2*b)/f1;
e=0; m=0; for(i=0;i<19;i++)
e=e+pow((a*exp((i+1)*b)-y[i]),2); for(i=0;i<18;i++){
y1[j]=log(y0[j]); x[j]=j+1; } y0 储存函数值，y1 储存变换后的值。 2、 此时可利用教材所给矩阵（4.61）进行计算，按照对应项构造系数矩阵，这一过程 程序代码为 m=19; p=q=s=t=0; for(k=0;k<19;k++){ p=p+x[k]; q=q+x[k]*x[k]; s=s+y1[k]; t=t+x[k]*y1[k]; } 3、 运用该矩阵解出我们所需要的 a1 与 b 的值，并变换 a 与 b 为原来的值，这一过程
（因偏导数方程过于繁琐，此处省略具体方程）。
接下来运用牛顿迭代法（弦割法）计算 a，b 的值。 （1）、先设定初始值 a0，b0； （2）、带入公式
f(a,b)=f(a0,b0)+f1(a0,b0)(a-a0)+f2(a0,b0)(b-b0)
g(a,b)=g(a0,b0)+g1(a0,b0)(a-a0)+g2(a0,b0)(b-b0)
可以看到，在选取合适的 a 与 b 的 值后，运用牛顿迭代法计算出的 a、b 估计值的精 确程度明显提高，可取数据 a=0.241187，b=0.308055（e=3.946625，m=5.234590）。 四、 结果说明 方法一中因存在计算机的舍入误差，对结果造成的影响较大；方法二的结果稍有改善， 不过因采取的牛顿迭代法，需要恰当的多次选取初值才有可能达到理想的目标。在我 看来，可以用第一种方法的值带入第二种方法进行适当的修正，如果达不到精度的要 求，再去寻找合适的初值。 五、 源程序 方法一、x1来自234
5
6
7
8
9
y 0.898 2.38 3.07 1.84 2.02 1.94 2.22 2.77 4,02
x 10 11 12 13 14 15 16 17 18 19
y 4.76 5.46 6.53 10.9 16.5 22.5 35.7 50.6 61.6 81.8
二、 方法简要（方法一为自主编制，方法二对网上的方法进行了参考）
printf("输入函数值"); for(i=0;i<19;i++){
scanf("%lf",&y0[i]); } for(j=0;j<19;j++){
y1[j]=log(y0[j]); x[j]=j+1; }
m=19; p=q=s=t=0; for(k=0;k<19;k++){
p=p+x[k]; q=q+x[k]*x[k]; s=s+y1[k]; t=t+x[k]*y1[k]; }
#include<stdio.h>
#include<math.h> double log(double x); main(void)
{ double x[19], y0[19], y1[19]; double p, q, m, s, t; double a0, a1, a, b, e=2.718282; double f, z, f0; int i, j, k, l;
if(pow((a*exp((i+1)*b)-y[i]),2)>m) m=pow((a*exp((i+1)*b)-y[i]),2);
} printf("a=%f b=%f 均方误差 e=%f 最大偏差 m=%f\n",a,b,e/19,sqrt(m)); } return 0;
}
printf("a=%lf,b=%lf\n 均方误差为%lf,最大偏差为%lf\n", a, b, f0, z); 注意，求均方误差与最大偏差要使用 y0 数组储存的函数值。 总体来说上述方法简便直观，但在计算数值时由于计算机存在舍入误差，可能会 对结果造成很大的影响。 5、 其次，我们可以不对方程进行变换，直接运用最小二乘法的定义，通过求解均方 误差方程的法方程组来确定式中 a、b 两个参数。这样会遇到一个很大的难题便是 非线性方程组的求解。参照网上给出的计算方法，我们采用弦割法（牛顿迭代法） 求解。 首先，设各结点的方差和为 e，有
for(i=0;i<19;i++) x[i]=i+1;