mathematica-数值命令与例题

第六章数值计算命令与例题北京交通大学6.1求近似函数•在生产和实验中,人们经常遇到需要通过某个未知的函数f(x)在有限个给定点的函数值:{xi,yi},i=1,2,….,n,这里f(xi)=yi去获得函数f(x)的近似函数(x),求近似函数(x)的方法主要有拟合方法和插值方法。6.1.1曲线拟合•曲线拟合主要用来求一元近似函数,它是根据最小二乘原理的意义下获得近似函数的,此近似函数具有在数据点处的误差平方和最小的特点。记函数集合:M=Span[0,1,2,…,m]={(x)|(x)=a00(x)+a11(x)+…+amm(x),aiR}•称集合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)=a00(x)+a11(x)+…+amm(x)形式的近似函数(x)命令形式2：Fit[{y1,y2,...,yn},{0,1,2,…,m},x]功能：根据数据点集{{1,y1},{2,y2},...,{n,yn}}求出具有拟合函数为(x)=a00(x)+a11(x)+…+amm(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)=a0+a1x+a2x2+…+amxm多项式拟合算法输入n+1个拟合点：(xi,yi),i=0,1,…,n根据散点图确定拟合多项式的次数m计算相应正规线性方程组的系数和右端项解正规正规线性方程组，得解：a0*,a1*,…,am*写出拟合多项式*(x)=a0*+a1*x+a2*x2+…+am*xm求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];说明：本程序用于求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}]是用简化的正规方程组编程的。例1．已知一组实验数据x1345678910f(x)1054211234用多项式拟合求其拟合曲线。解：执行m次多项式拟合程序后，在输入的两个窗口中按提示分别输入{1,3,4,5,6,7,8,9,10}，{10,5,4,2,1,1,2,3,4}每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上画出散点图。由于该散点图具有2次多项式形状，因此在确定选择多项式次数窗口输入2，按OK”按扭后得如下输出结果。a=953381533813017381301725317b=32147102513.4597-3.60531t+0.267571t2于是得求出的二次多项式拟合函数为(t)=13.4597-3.60531t+0.267571t2而且从图形上看拟合效果很好。例2．已知一组实验数据x681012141618202224f(x)4.64.84.64.95.05.45.15.55.66.0修改多项式拟合程序使其具有可以选择不同多项式拟合函数功能，并用此程序确定本题多项式拟合曲线。解：在拟合多项式程序中加入评价拟合效果的判定人机交互语句tu=Input[Satisfatory?0(No)or1{Yes}]和While语句来调控何时终止实验，调控变量用tu取值确定，取值0表示不满意，1表示满意。此外，去掉正规方程组系数及拟合多项式的输出，代之以在图形上表出拟合多项式的次数提示，修改后的程序如下。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拟合多项式},DisplayFunction-Identity];Show[{p1,h},DisplayFunction-$DisplayFunction];tu=Input[Satisfatory?0(No)or1{Yes}]]执行该程序后，屏幕上出现拟合多项式次数输入窗口和散点图。由于该散点图不好确定拟合次数，先用3次拟合多项式计算，因此输入“3”并用鼠标点击窗口的“OK”按扭，得如下输出图形。屏幕出现提示是否满意的输入窗口，因为不太满意，输入：0，单击“OK”按扭并在屏幕上出现拟合多项式次数输入窗口中输入：6，单击OK，屏幕出现下一个组合图形。继续做实验，得到如下若干图形。由于总共有10个数据点，所以拟合多项式最高次数只能到9次，因此实验到9次拟合多项式后，在屏幕出现提示是否满意的输入窗口中，输入：1，单击“OK”按扭退出实验。线性模型拟合线性模型拟合算法1.输入n+1个拟合点：(xi,yi),i=0,1,…,n2.根据散点图确定拟合基函数组3.计算相应正规线性方程组的系数和右端项4.解正规正规线性方程组，得解：a0*,a1*,…,am*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]说明：本程序用于求线性模型拟合。程序执行后，按要求通过键盘输入插值基点xi:{x0,x1,...,xn}、对应函数值yi:{y0,y1,…,yn}后，计算机给出散点图和请求输入拟合拟合基函数组{0(x)，1(x)，2(x)、…、m(x)}的窗口，操作者可以根据散点图确定拟合基函数组通过键盘输入，程序即可给出对应的正规方程组系数矩阵a、常数项b、线性模型拟合函数和由拟合函数图形与散点图画在一起的图形。程序中变量说明：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)=a00(x)+a11(x)+…+amm(x)形式的近似函数(x)。例3．已知函数在自变量x=1,2,…,10上数据为{2.89229,2.86323,0.473147,-2.25209,-2.87003,-