第六章 数值计算命令与例题

6.1.1 曲线拟合 • 曲线拟合主要用来求一元近似函数, 它是根据最小二乘原 理的意义下获得近似函数的, 此近似函数具有在数据点处的 误差平方和最小的特点。记函数集合:
M=Span[0, 1, 2,…, m]= {(x)|(x)= a0 0(x)+a11(x)+…+amm(x), ai R}

xi=Table[i,{i,6,24,2}]; yi={4.6,4.8,4.6,4.9,5,5.4,5.1,5.5,5.6,6}; n=Length[xi]; h=ListPlot[Table[{xi[[i]],yi[[i]]},{i,1,n}],PlotStyle->PointSize[0.04]] tu=0; While[tu==0, m=Input["多项式次数m"]; s=Table[Sum[xi[[k]]^i,{k,1,n}],{i,0,2m}]; a=Table[s[[i+j-1]],{i,1,m+1},{j,1,m+1}]; (*Print["a=",MatrixForm[a]];*) b=Table[Sum[xi[[k]]^i*yi[[k]],{k,1,n}],{i,0,m}]; (*Print["b=",MatrixForm[b]];*) xx=Table[x[i],{i,1,m+1}]; g=Solve[a.xx==b,xx]; fa=Sum[x[i]*t^(i-1),{i,1,m+1}]/.g[[1]]; p=fa//N; p1=Plot[p,{t,xi[[1]],xi[[n]]},PlotLabel->{m"拟合多项式"},
•
称集合M为函数0, 1, 2,…, m张成的空间，m+1个函 数0(x), 1(x), 2(x),…, m(x)称为拟合基函数集合, 它们都是已知的函数。
Mathematica曲线拟合的一般形式为: Fit[{数据点集合}, {拟合基函数集合}, 自变量名] 具体的拟合命令有: 命令形式1：Fit[{{x1,y1},{x2,y2},...,{xn,yn}},{ 0, 1, 2,…, m },x] 功能：根据数据点集{{x1,y1},{x2,y2},...,{xn,yn}}求出具有拟合函数为 (x)= a 0 0(x)+a11(x)+…+a mm(x) 形式的近似函数(x)

说明：本程序用于求m次多项式拟合。程序执 行后，按要求通过键盘输入拟合基点xi:{x0 , x1, ... , xn }、对应函数值yi:{ y0 , y1 , … , yn } 后，计算机给出散点图和请求输入拟合多项式 次数的窗口，操作者可以根据散点图确定拟合 多项式的次数通过键盘输入，程序即可给出对 应的正规方程组系数矩阵a、常数项b、m次拟 合多项式和由拟合函数图形和散点图画在一起 的图形。






程序中变量说明 xi：存放拟合基点{x0 , x1, ... , xn } yi: 存放对应函数值{y0 , y1 , … , yn} m: 存放拟合多项式次数 a: 存放正规方程组系数矩阵 b: 存放正规方程组常数项 p: 存放m次拟合多项式 h: 存放散点图 p1：存放拟合函数图形 xx:定义正规方程组变量，存放m次拟合多项式的系数 注：语句s=Table[Sum[xi[[k]]^i,{k,1,n}],{i,0,2m}]、 a=Table[s[[i+j-1]],{i,1,m+1},{j,1,m+1}]、 b=Table[Sum[xi[[k]]^i*yi[[k]],{k,1,n}],{i,0,m}]是用简化 的正规方程组编程的。
求m次多项式拟合程序


Clear[xi,xx,yi]; xi=Input["xi="] yi=Input["yi="] n=Length[xi]; h=ListPlot[Table[{xi[[i]],yi[[i]]},{i,1,n}],PlotStyle->PointSize[0.04]] m=Input["多项式次数m="] s=Table[Sum[xi[[k]]^i,{k,1,n}],{i,0,2m}]; a=Table[s[[i+j-1]],{i,1,m+1},{j,1,m+1}]; Print["a=",MatrixForm[a]]; b=Table[Sum[xi[[k]]^i*yi[[k]],{k,1,n}],{i,0,m}]; Print["b=",MatrixForm[b]]; xx=Table[x[i],{i,1,m+1}]; g=Solve[a.xx==b,xx]; fa=Sum[x[i]*t^(i-1),{i,1,m+1}]/.g[[1]]; p=fa//N p1=Plot[p,{t,xi[[1]],xi[[n]]},DisplayFunction->Identity]; Show[{p1,h},DisplayFunction->$DisplayFunction];
例1．已知一组实验数据 x 1 3 4 5 6 7 8 9 10 f(x) 10 5 4 2 1 1 2 3 4 用多项式拟合求其拟合曲线。
解：执行m次多项式拟合程序后，在输入的两 个窗口中按提示分别输入 {1,3,4,5,6,7,8,9,10}，{10,5,4,2,1,1,2,3,4} 每次输入后用鼠标点击窗口的“OK”按扭，计 算机在屏幕上画出散点图。

由于该散点图不好确定拟合次数，先用3次拟合多项式计 算，因此输入“3”并用鼠标点击窗口的“OK”按扭，得如 下输出图形。
屏幕出现提示是否满意的输入窗口，因为不太满意，输入： 0，单击“OK”按扭并在屏幕上出现拟合多项式次数输入 窗口中输入：6，单击OK，屏幕出现下一个组合图形。
继续做实验，得到如下若干图形。

程序中变量说明： xi：存放拟合基点{x0 , x1, ... , xn } yi: 存放对应函数值{y0 , y1 , … , yn} m1: 存放拟合基函数组{0(x)，1(x)，2(x)、…、m(x)} m: 存放拟合基函数组个数 a: 存放正规方程组系数矩阵 b: 存放正规方程组常数项 p: 存放拟合基函数组在拟合基点{x0 , x1, ... , xn }的函数值 pp: 存放求出的线性模型拟合函数 h: 存放散点图 p1：存放拟合函数图形 xx:定义正规方程组变量，存放线性模型拟合系数 注：（1）语句m1=Input["输入拟合基函数组"]产生的输入应用函数表，即用 “{0[x]，1[x]，…，m[x]}”的形式输入。 （2）Mathematica数学软件中也有一个求线性模型拟合的命令，命令形式为 Fit[{{x1,y1},{x2,y2},...,{xn,yn}},{ 0, 1, 2,…, m },x] 它可以根据数据点集{{x1,y1},{x2,y2},...,{xn,yn}}求出具有拟合函数为 (x)= a 0 0(x)+a11(x)+…+a mm(x) 形式的近似函数(x)。
说明：本程序用于求线性模型拟合。程序 执行后，按要求通过键盘输入插值基点 xi:{x0 , x1, ... , xn }、对应函数值yi:{ y0 , y1 , … , yn }后，计算机给出散点图和请求 输入拟合拟合基函数组{0(x)，1(x)， 2(x)、…、m(x)}的窗口，操作者可以根 据散点图确定拟合基函数组通过键盘输入， 程序即可给出对应的正规方程组系数矩阵a、 常数项b、线性模型拟合函数和由拟合函数 图形与散点图画在一起的图形。
命令形式2：Fit[{y1,y2,...,yn},{ 0, 1, 2,…, m },x] 功能：根据数据点集{{1,y1},{2,y2},...,{n,yn}}求出具有拟合函数为 (x)= a 0 0(x)+a11(x)+…+a mm(x) 形式的近似函数(x) 命令形式3：Fit[{{x1,y1},{x2,y2},...,{xn, yn}}, Table[x^i,{i,0,m}] ,x] 功能：根据数据点集{{x1,y1},{x2,y2},...,{xn,yn}}求出拟合函数为m次多 项式的近似函数(x) =a 0 +a1x+ a2x2 +…+a mx m
由于总共有10个数据点，所以拟合多项式最高次 数只能到9次，因此实验到9次拟合多项式后，在 屏幕出现提示是否满意的输入窗口中，输入：1， 单击“OK”按扭退出实验。
线性模型拟合
线性模型拟合算法 1.输入n+1个拟合点： (xi, yi),i=0,1,…,n 2.根据散点图确定拟合基函数组 3.计算相应正规线性方程组的系数和右端项 4.解正规正规线性方程组，得解：a0*,a1*,…,a m* 5.写出线性拟合模型*(x)= a0*0(x)+ a1*1(x)+ a2*2(x)+ …+ am*m(x)
求线性模型拟合程序
Clear[x,xi,xx,yi]; xi=Input["xi="] yi=Input["yi="] n=Length[xi]; h=ListPlot[Table[{xi[[i]],yi[[i]]},{i,1,n}],PlotStyle->PointSize[0.04]] m1=Input["输入拟合基函数组"] m=Length[m1] p=Table[m1/.x->xi[[k]],{k,1,n}] a=Table[Sum[p[[k,i]]*p[[k,j]],{k,1,n}],{i,1,m},{j,1,m}]//N (*Print["a=",MatrixForm[a]];*) b=Table[Sum[p[[k,i]]*yi[[k]],{k,1,n}],{i,1,m}]//N (*Print["b=",MatrixForm[b]];*) xx=Table[xt[i],{i,1,m}] g=Solve[a.xx==b,xx] fa=Sum[xt[i]*m1[[i]],{i,1,m}]/.g[[1]] pp=fa//N; p1=Plot[pp,{x,xi[[1]],xi[[n]]},DisplayFunction->Identity]; Show[{p1,h},DisplayFunction->$DisplayFunction]