工数上机.doc

一、牛顿迭代法（C语言） （1）介绍 设r是f(x) = 0的根，选取x0作为r初始近似值，过点（x0,f(x0)）做曲线y = f(x)的切线L，L的方程为y = f(x0)+f(x0)(x-x0)，求出L与x轴交点的横坐标 x1 = x0-f(x0)/f(x0)，称x1为r的一次近似值。过点（x1,f(x1)）做曲线y = f(x)的切线，并求该切线与x轴交点的横坐标 x2 = x1-f(x1)/f(x1)，称x2为r的二次近似值。重复以上过程，得r的近似值序列，其中x(n+1)=x(n)－f(x(n))/f(x(n))，称为r的n+1次近似值，上式称为牛顿迭代公式。 　　 解非线性方程f(x)=0的牛顿法是把非线性方程线性化的一种近似方法。把f(x)在x0点附近展开成泰勒级数 f(x) = f(x0)+(x－x0)f(x0)+(x－x0)^2*f(x0)/2! +… 取其线性部分，作为非线性方程f(x) = 0的近似方程，即泰勒展开的前两项，则有f(x0)+f(x0)(x－x0)=f(x)=0 设f(x0)≠0则其解为x1=x0－f(x0)/f(x0) 这样，得到牛顿法的一个迭代序列：x(n+1)=x(n)－f(x(n))/f(x(n))。 （2）程序流程分析： 赋值x0＝1.5,即迭代初值； 用初值x0代入方程中计算此时的f(x0)及f’(x0)，程序中用变量f描述方程的值，用fd描述方程求导之后的值； 计算增量d=f/fd； 计算下一个x，x=x0-d; 把新产生的x替换x0，为下一次迭代做好准备; ⑥ 若d绝对值大于1e-，则重复②③④⑤步。#include stdio.h #include math.h typedef double (*Fun)(double); double newtown(double x0, Fun f, Fun df) { double x1 = x0 - f(x0) / df(x0); while(fabs(x0 - x1) 1e-2) { x0 = x1; x1 = x0 - f(x0) / df(x0); } return x1; } double fun(double x) { double x2 = x * x, x3 = x2 * x; return x3-2*x-5; } double dfun(double x) { double x2 = x * x * x; return 3*x2 - 2; } int main(void) { /* int i = 0; double x = 0.0; for( ; i 20; i++) { printf(f(%.10f) = %.10f\n, x, fun(x)); x += 0.1; }*/ double ret = newtown(1.5, fun, dfun); printf(%.10f\n, ret); printf(f(%.10f) = %.10f\n, ret, fun(ret)); return 0; } （4）结果： 二、最小二乘法拟合曲线（matlab） （1）程序：由于此处只需要一个X对应一个Y,则不需要输入X2和X3的值 x1=[0.5 1.5 2.5 3.5 4.5 5.5]; x2=[]; x3=[]; y=[0.0001 1.9998 4.0010 5.9980 8.0001 10.00010]; format short g Y=y X11=[ones(1,length(y));x1;x2;x3] B1=regress(Y,X11)% 多元一次线性回归 [m,n]=size(X11) X22=[]; for i=2:n for j=2:n if i=j X22=([X22,X11(:,i).*X11(:,j)]); else continue end end end X=[X11,X22]; B2=regress(Y,X)% 多元二次线性回归 [Y X*B2 Y-X*B2] plot(Y,X11*B1,o,Y,X*B2,*) hold on,line([min(y),max(y)],[min(y),max(y)]) axis([min(y) max(y) min(y) max(y)]) legend(一次线性回归,二次线性回归) xlabel(实际值);ylabel(计算值) （2）运算结果： Y = 0.0001 1.99