第六章 数值计算命令与例题 北京交通大学

{(x)|(x)= a0 0(x)+a11(x)+…+amm(x), ai R} 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)= a 0 0(x)+a11(x)+…+a mm(x) 形式的近似函数(x) 命令形式2：Fit[{y1,y2,...,yn},{ 0, 1, 2,…, m },x] 功能：根据数据点集{{1,y1},{2,y2},...,{n,yn}}求出具有拟合函数为 命令形式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

多项式拟合算法 输入n+1个拟合点： (xi, yi),i=0,1,…,n 根据散点图确定拟合多项式的次数m 计算相应正规线性方程组的系数和右端项 解正规正规线性方程组，得解：a0*,a1*,…,a m* 写出拟合多项式*(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, 说明：本程序用于求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．已知一组实验数据 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”按扭，计算机在屏幕上画出散点图。

由于该散点图具有2次多项式形状，因此在确定选择多项式次数窗口输入2，按OK”按扭后得如下输出结果。 a= 9 53 381 53 381 3017 381 3017 25317 b=32 147 1025 13.4597 - 3.60531 t + 0.267571 t2

于是得求出的二次多项式拟合函数为 (t)=13.4597 - 3.60531 t + 0.267571 t2 而且从图形上看拟合效果很好。

例 2．已知一组实验数据 x 6 8 10 12 14 16 18 20 22 24 f(x) 4.6 4.8 4.6 4.9 5.0 5.4 5.1 5.5 5.6 6.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*,…,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]

说明：本程序用于求线性模型拟合。程序执行后，按要求通过键盘输入插值基点xi:{x0 , x1, 说明：本程序用于求线性模型拟合。程序执行后，按要求通过键盘输入插值基点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)= a 0 0(x)+a11(x)+…+a mm(x) 形式的近似函数(x)。

例3．已知函数在自变量x=1,2, …, 10上数据为{2.89229, 2.86323, 0.473147, -2.25209, -2.87003, -0.835768, 1.97187, 2.96841, 1.23648, -1.63202}，试用合适的函数进行拟合。 解：执行线性模型拟合程序后，在输入的两个窗口中按提示分别输入 {1,2,3,4,5,6,7,8,9,10}、{2.89229, 2.86323, 0.473147, -2.25209, -2.87003, -0.835768, 1.97187, 2.96841, 1.23648, -1.63202} 每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上画出散点图。

由于该散点图具有类似正弦曲线的形状，因此在确定选择拟合基函数窗口输入“{1,Sin[x]}”，按“OK”按扭后得如下输出结果。 a=10. 1.41119 1.41119 5.00143 b=4.81552 15.4239 0.0482789 + 3.07027 Sin[x]

于是，我们得到了很好的拟合函数(x)=0.0482789 + 3.07027 sin(x) 。 对于本题，如果看到散点图具有一个弯曲而选用三次多项式拟合，则输入“{1,x,x^2,x^3}”后会得到如下输出。

a= 10. 55. 385. 3025. 55. 385. 3025. 25333. 385. 3025. 25333. 220825. 3025. 25333. 220825. 1.97841 106 b=4.81552 14.0236 104.284 820.711 10.4765 - 7.37686 x + 1.41989 x2 - 0.0796299 x3

显然这个拟合很不满意。

6.1.2 函数插值 多项式插值 多项式插值是在给定n个数据点:{xi, yi}, i=1,2,…., n 后, 求出一个次数不超过n-1的多项式(x)作为函数 y=f(x) 的近似表达式，它就是计算方法中常说的拉格朗日(Lagrange)插值或Newton插值多项式。 Mathematica多项式插值的一般形式为: InterpolatingPolynomial [{数据点集合}, 自变量名] 具体的多项式插值命令有: 命令形式1：InterpolatingPolynomial [{{x1,y1},{x2,y2},...,{xn,yn}},x] 功能：根据数据点集{{x1,y1},{x2,y2},...,{xn,yn}}求出n-1次插值多项式(x) 命令形式2：InterpolatingPolynomial [{y1,y2,...,yn},x] 功能：根据数据点集{{1,y1},{2,y2},...,{n,yn}}求出n-1次插值多项式(x)

Lagrange插值 求n次Lagrange插值多项式算法 输入n+1个插值点： (xi, yi),i=0,1,…,n 计算插值基函数l0n(x), l1n(x) ,…, lnn(x) 给出n次Lagrange插值多项式：Ln(x)= y0 l0n(x)+ y1 l1n(x) +…+yn lnn(x)

求Lagrange插值多项式程序 Clear[lag,xi,x,yi]; xi=Input["xi="] yi=Input["yi="] n=Length[xi]-1; p=Sum[yi[[i]]*(Product[(x-xi[[j]])/(xi[[i]]-xi[[j]]),{j,1,i-1}] *Product[(x-xi[[j]])/(xi[[i]]-xi[[j]]),{j,i+1,n+1}]),{i,1,n+1}]; lag[x_]=Simplify[p] 说明：本程序用于求n次Lagrange插值多项式。程序执行后，按要求通过键盘输入插值基点xi:{x0 , x1, ... , xn }和对应函数值yi:{ y0 , y1 , … , yn }后，程序即可给出对应的n次Lagrange插值多项式lag[x]。

程序中变量说明 xi：存放插值基点{x0 , x1, ... , xn } yi: 存放对应函数值{y0 , y1 , … , yn} lag[x]: 存放求出的n次Lagrange插值多项式Ln(x) 注：语句lag[x_]=Simplify[p]用简化形式给出对应的n次Lagrange插值多项式。

例．给定数据表 x 0 1 2 3 y=f(x) 1 3 5 12 用Lagrange插值法求三次插值多项式,并给出函数f(x)在x =1.4的近似值。

解： 执行Lagrange插值程序后，在输入的两个窗口中按提示分别输入{0, 1, 2, 3}、{1, 3, 5, 12}，每次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数。 6 + 22 x - 15 x2 + 5 x3 ----------------------- 6 所以得到三次插值多项式L3(x)=1+11 x/3-5 x2/2+5 x3/6 接着键入“lag[1.4]”，则输出3.52，因此f(x)在x =1.4的近似值为3.52，即f(1.4)3.52.

Newton插值 求Newton插值多项式算法 1. 输入n+1个插值点： (xi, yi),i=0,1,…,n 2. 计算差商表

求Newton插值多项式程序 Clear[newt,s,x]; xi=Input["xi="] yi=Input["yi="] n=Length[xi]; (*计算差商表*) f=Table[0,{n},{n}]; Do[f[[i,1]]=yi[[i]],{i,1,n}] Do[f[[i,j+1]]=(f[[i,j]]-f[[i+1,j]])/(xi[[i]]-xi[[i+j]]), {j,1,n-1},{i,1,n-j}] Print["差商表"] Do[Print[xi[[i]]," ",f[[i]]],{i,1,n}] (*求Newton插值多项式*) fa=1; s=f[[1,1]]; Do[fa=(x-xi[[k]])*fa;s=s+fa*f[[1,k+1]],{k,1,n-1}] newt[x_]=s Simplify[%]

说明：本程序用于求n次Newton插值多项式。程序执行后，按要求通过键盘输入插值基点xi:{x0 , x1, 说明：本程序用于求n次Newton插值多项式。程序执行后，按要求通过键盘输入插值基点xi:{x0 , x1, ... , xn }和对应函数值yi:{ y0 , y1 , … , yn }后，程序依次给出输入的数据表、计算出的差商表、Newton插值多项式、Newton插值多项式的简化形式。 程序中变量说明： xi：存放插值基点{x0 , x1, ... , xn } yi: 存放对应函数值{y0 , y1 , … , yn} f:存放函数值{y0 , y1 , … , yn}及所有差商 newt[x]: 存放求出的n次newton插值多项式Nn(x) 注： （1)语句f=Table[0,{n},{n}]用于产生一个nn的矩阵变元用于存放函数值{y0 , y1 , … , yn}及所有差商。 （2)在Mathematica中有一个求n次插值多项式的命令，命令形式 InterpolatingPolynomial[{{x01,y0}，{x1,y1},{x2,y2},…,{xn,yn}}, x] 它可以求过n+1个插值点{{x01,y0}，{x1,y1},{x2,y2},…,{xn,yn}}的n次插值多项式Pn(x)。

例2．给定数据表 x 4.0002 4.0104 4.0233 4.0294 y 0.6020817 0.6031877 0.6045824 0.6052404 (1)计算差商表 (2)用Newton插值法求三次插值多项式Nn(x) (3)求f(4.011)的近似值 解： 执行Newton插值程序后，在输入的两个窗口中按提示分别输入 {4.0002, 4.0104, 4.0233, 4.0294}，{0.6020817, 0.6031877, 0.6045824, 0.6052404}, 每次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数。

{4.0002, 4.0104, 4.0233, 4.0294} {0.6020817, 0.6031877, 0.6045824, 0.6052404} 差商表 4.0002 {0.6020817, 0.108431, -0.0136404, 0.0211629} 4.0104 {0.6031877, 0.108116, -0.0130225, 0} 4.0233 {0.6045824, 0.107869, 0, 0} 4.0294 {0.6052404, 0, 0, 0} 0.6020817 + 0.108431 (-4.0002 + x) - 0.0136404 (-4.0104 + x) (-4.0002 + x) + 0.0211629 (-4.0233 + x) (-4.0104 + x) (-4.0002 + x) -1.41642 + 1.23926 x - 0.268313 x2 + 0.0211629 x3

于是我们得到本题的差商表为 xi yi f[, ] f[ , , ] f[ , , ,] 4.0002 0.6020817, 0.108431, -0.0136404, 0.0211629 4.0104 0.6031877, 0.108116, -0.0130225 4.0233 0.6045824, 0.107869 4.0294 0.6052404 和三次插值多项式 N3(x)= -1.41642 + 1.23926 x - 0.268313 x2 + 0.0211629 x3 接着键入“newt[4.011]”，则输出0.603253，因此f(4.011) 0.603253。

例3．多项式插值的误差估计式中可以看到，当插值节点越多时误差会越小，这个结论正确吗？通过实验说明该结论的正确性。 解: 考虑函数f(x) = (1+x2)-1 在区间[-4,4]内选取不同个数的等距插值节点做观察，这里分别选[-4,4] 内的9个和11个的等距节点来做实验，将对应的插值函数图与被插函数f(x) = (1+x2)-1画在一起做观察，为简单起见，这里用Mathematica 命令做实验，对应命令为

In[1]:= u=Table[{x,(1+x^2)^-1},{x,-4,4}] ; （*采取f(x) 在[-4,4] 内的9个插值点 In[2]:= g=ListPlot[u, PlotStyle->PointSize[0.04]] （*将散点图图形文件存放在变量g中 In[3]:= s=InterpolatingPolynomial[u , x] ; （*将插值函数存放在变量s中 In[4]:= t= Plot[{s, (1+x^2)^-1}, {x,-4,4}, PlotStyle->{{Thickness[0.005]},{Thickness[0.006]}}] （*将插值函数s与f(x)画在一起的图形文件存放在变量t中 In[5]:= Show[t, g]（*将散点图, 插值函数s, f(x)画在一个坐标系中

在[-4,4]中选9个等距节点的插值函数与被插函数图，粗线为被插函数图

In[6]:= u=Table[{x,(1+x^2)^-1},{x,-4,4，0. 8}] ; （ In[6]:= u=Table[{x,(1+x^2)^-1},{x,-4,4，0.8}] ; （*采取f(x) 在[-4,4] 内的11插值点 In[7]:= g=ListPlot[u, PlotStyle->PointSize[0.04]] In[8]:= s=InterpolatingPolynomial[u , x] ; In[9]:= t= Plot[{s, (1+x^2)^-1}, {x,-4,4}, PlotStyle->{{Thickness[0.005]},{Thickness[0.006]}}] In[10]:= Show[t, g]

在[-4,4]中选11个等距节点的插值函数与被插函数图

Interpolation [{数据点集合}] 分段多项式插值 Mathematica分段多项式插值的一般形式为: Interpolation [{数据点集合}] 具体的分段多项式插值命令有: 命令形式1：Interpolation [{{x1,y1},{x2,y2},...,{xn,yn}}] 功能：根据数据点集{{x1,y1},{x2,y2},...,{xn,yn}}求出分段插值多项式(x) 命令形式2：Interpolation [{y1,y2,...,yn}] 功能：根据数据点集{{1,y1},{2,y2},...,{n,yn}}求出分段插值多项式(x)

例5: 用分段多项式插值重做例4的插值问题, 并验证, 随着插值点的增多, 分段插值函数的误差会越来越小。(见图) 解: In[27]:= r=Interpolation[Table[{x,(1+x^2)^-1}, {x,-4,4}] ] Out[27]= InterpolatingFunction[{-4, 4}, <>] In[28]:= (1+x^2)^-1/.x->3.7 Out[28]= 0.0680735 In[29]:= r[3.7] Out[29]=0.0734 In[30]:= d= Table[{x,(1+x^2)^-1}, {x,-4,4,0.5}]; In[31] := r1=Interpolation[d] Out[31]= InterpolatingFunction[{-4, 4.}, <>] In[32]:= r1[3.7] Out[32]= 0.0681761 In[33]:= Plot[{r1[x], (1+x^2)^-1}, {x,-4,4}] 下一页

返回

分段线性插值 用分段线性插值函数做近似计算的算法 1. 输入n个插值点： (xi, yi),i=1,…,n 2. 输入要近似计算的自变量点xa 3.寻找包含xa的小区间[xi , xi+1] 4.用 [xi , xi+1] 上的线性插值函数在xa处的值作为f(xa)的函数值

用分段线性插值函数做近似计算的程序 Clear[x,a,b]; li[a_,b_,x_]:=(b[[2]]-a[[2]])/(b[[1]]-a[[1]])*(x-a[[1]])+a[[2]] xi=Input["xi="] yi=Input["yi="] xa=Input["xa="] n=Length[xi]; If[xa<xi[[1]]||xa>xi[[n]],Print["超限"];Break[], Do[If[xa>=xi[[k]],m=k],{k,1,n-1}]; Print["m=",m,"x=",xi[[m]]]; li[{xi[[m]],yi[[m]]},{xi[[m+1]],yi[[m+1]]},xa]] 说明：本程序用分段线性插值做近似计算。程序执行后，按要求通过键盘输入插值基点xi:{ x1, ... , xn }和对应函数值yi:{ y1 , … , yn }和要做近似计算的点xa后，程序将给出包含xa小区间[xi , xi+1]的下标i和xi ,和函数f(xa)的近似值。如果xa不在[x1, xn ]内程序给出超限提示。

程序中变量说明 xi：存放插值基点{ x1, ... , xn } yi: 存放对应函数值{ y1 , … , yn} xa: 存放要做近似计算的自变量值 m:存放包含xa小区间[xi , xi+1]的下标i 注：在Mathematica中有一个求分段线性插值函数的命令，命令形式为 Interpolation [{{x1,y1},{x2,y2},...,{xn,yn}}，InterpolationOrder -> 1] 用如上Mathematica命令求出的分段线性插值函数没有给出具体的分段函数表达式, 而是用 “InterpolatingFunction[{x1, xn}, <>] ”作为所求的分段插值函数, 通常可以用 变量=Interpolation [{数据点集合}] 把所得的分段插值函数存放在变量中, 如果要计算插值函数在某一点的如p点的值,只要输入 “变量[p]” 即可，如果要表示这个插值函数应该用 “变量[x]” 。 此外，命令 Interpolation [{{x1,y1},{x2,y2},...,{xn,yn}}，InterpolationOrder -> 3] 可以给出更好的分段插值函数。

例4．给定数据表 xi 1 2 3 4 yi 0 -5 -6 3 （1）用分段线性插值函数求函数f(x)在x =1.4的近似值 （2）用Mathematica命令画出分段线性插值图形

解：执行求分段线性插值函数程序后，在输入的两个窗口中按提示分别输入 {1，2，3，4}，{0，-5，-6，3}，1.4 每次输入后用鼠标点击窗口的“OK”按扭，得如下插值函数。 m=1 x=1 -2 此结果说明f(x)在x =1.4的近似值为2，它对应的小区间为[x1 , x2]。为完成第二个问题，键入命令 d=Interpolation [{{1,0},{2,-5},{3,-6},{4,3}},InterpolationOrder -> 1] Plot[d[x],{x,1,4}]

得输出如下。 InterpolatingFunction[{1, 4}, <>]

6.2 解常微分方程 自变量、未知函数以及未知函数的导数（或微分）之间的一个（或一组）关系式称为微分方程（或微分方程组）。求出常微分方程（或微分方程组）的未知函数（或未知函数组），称为解常微分方程。含有任意常数的解称为通解，否则称为特解。n阶常微分方程的一般形式为： 利用Mathematica可以象高等数学一样求出常微分方程的解析解（准确解），如果求不出解析解Mathematica可以求出微分方程的数值解（即近似解。Mathematica解常微分方程的命令有求常微分方程（组）的解析解命令和求常微分方程（组）的数值解命令。

6.2.1求常微分方程（组）的解析解 命令形式1：DSolve[eqn, y[x], x] 命令形式2：DSolve[{eqn1,eqn2, ... }, {y1[x],y2[x], ... }, x] 功能：求出常微分方程组{eqn1,eqn2, ... }的所有未知函数{y1[x],y2[x], ... }的解析通解。 注意： 每个常微分方程的中的等号输入时应该用两个等号代替 未知函数及未知函数的导数要用“[自变量名]”指示未知函数的自变量 常微分方程组和未知函数组应该用大括号括起来 命令形式2可以用来求常微分方程的特解 eqn表示常微分方程, x是自变量名, y[x] 表示未知函数, 自变量名和未知函数可以是其他的变量名

例6: 求常微分方程y =2ax的通解，a为常数 解:Mathematica 命令: In[34]:= DSolve[y'[x] == 2 a x, y[x], x] Out[34]= {{y[x] -> a x 2 + C[1]}} 即得本问题的通解 例7: 求常微分方程y +y =0的通解 In[35]:= DSolve[y''[x]+y[x]==0, y[x], x] Out[35]= {{y[x] -> C[2] Cos[x] - C[1] Sin[x]}} 即得本问题的通解: C1，C2是任意常数。

例8: 求常微分方程 的特解。 解:Mathematica 命令: In[36]:= DSolve[{y'[x]==x/y[x]+y[x]/x, y[1]==2}, y[x], x] Out[36]= {{y[x] -> Sqrt[x2 (4 + 2 Log[x])]}} 本问题的特解为: 下面的操作可以检验所求特解的正确性: In[37]:= y[x_]:= Sqrt[x^2* (4 + 2 Log[x])] In[38]:= y[1] Out[38]=2 In[39]:=D[y[x],x]-x/y[x]-y[x]/x; In[40]:= Simplify[%] Out[40]= 0

命令形式1：NDSolve[ eqns, y, {x, xmin, xmax}] 例9: 求常微分方程组y(t) =x(t), x(t) =y(t) 的通解。 解:Mathematica 命令: In[41]:= DSolve[{y'[t]==x[t],x'[t]==y[t]}, {x[t],y[t]},t] C[1]+ E 2 t C[1]-C[2]+ E 2 t C[2] Out[41]= {{x[t] -> ------------------------------------------, 2 E t -C[1]+ E 2 t C[1]+C[2]+ E 2 t C[2] y[t] -> ----------------------------------------------}} 6.2.2 求常微分方程（组）的数值解命令 命令形式1：NDSolve[ eqns, y, {x, xmin, xmax}] 功能：求出自变量范围为[xmin, xmax]且满足给定常微分方程及初值条件eqns的未知函数y的数值解 命令形式2：NDSolve[ eqns, {y1,y2, ... }, {x, xmin, xmax}] 功能：求出自变量范围为[xmin, xmax]且满足给定常微分方程及初值条件eqns的未知函数y1,y2, ... 的数值解。

例10:求微分方程初值问题y=x+y,y(0)=1在区间[0,0.5]内的数值解 解:Mathematica 命令: In[42]:= NDSolve[ {y'[x]==x+y[x], y[0]==1}, y, {x,0,0.5} ] Out[42]= {{y -> InterpolatingFunction[{0., 0.5}, <>]}} 上面显示的是所求数值解的替换形式，为得到本问题数值解，再键入: In[43]:= f= y/.% [ [1] ] Out[43]= InterpolatingFunction[{0., 0.5}, <>] In[44]:= Plot[f[x],{x,0,0.5}] In[45]:= Table[ f[x], {x, 0, 0.5, 0.1} ] Out[45]= {1., 1.11034, 1.24281, 1.39972, 1.58365, 1.79744}

求函数x(t)和y(t)在[0,10]上的数值解。 例11: 已知常微分方程组 求函数x(t)和y(t)在[0,10]上的数值解。 解:Mathematica 命令: In[46]:=q=NDSolve[{x'[t]==-x[t]^2-y[t],y'[t]==2x[t]-y[t],x[0]==y[0]==1},{x,y},{t,0,10}] Out[46]= {{x -> InterpolatingFunction[{0., 10.}, <>], y -> InterpolatingFunction[{0., 10.}, <>]}} In[47]:= f=x/. q[[1]]; g=y/. q [[1]]; In[48]= f[0.4] Out[48]= 0.375078 In[49]= g[4] Out[49]= -0.0912007 In[50]= ParametricPlot[{f[t],g[t]},{t,0,10}]

Euler方法 Euler方法算法 输入函数f(x,y)、初值y0、自变量区间端点a,b步长h 计算节点数n和节点x k 用Euler公式y n+1= yn+h f(xn,yn) 求数值解

Euler方法程序 Clear[x,y,f] f[x_,y_]= Input["函数f(x,y)="] y0=Input["初值y0 ="] a=Input["左端点a="] b=Input["右端点b="] h=Input["步长h="] n=(b-a)/h; For[i=0,i<n,i++, xk=a+i*h; y1=y0+h*f[xk,y0]; Print["y(",xk+h//N,")=",y1//N]; y0=y1] 说明：本程序用Euler公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。

程序中变量说明 f[x,y]: 存放函数f(x,y) y0: 存放初值y0及数值解 a：存放自变量区间左端点 b: 存放自变量区间右端点 n: 存放节点个数 h: 存放节点步长 xk：存放节点xi y1: 存放数值解 注：（1）语句Print["y(",xk+h//N,")=",y1//N]是将数值解用6位有效数字显示出来，如果要显示n位数有效数字，可以将语句改为Print["y(",xk+h//N,")=",N[y1,n]] （2）Mathematica中有求微分方程初值问题数值解的命令，形式为: NDSolve[ {y’[x]==f[x,y[x]],y[a]==y0}, y, {x, a, b}] 由NDSolve命令得到的解是以{{未知函数名->InterpolatingFunction[range, <>]}}的形式给出的,其中的InterpolatingFunction[range, <>]是所求的插值函数表示的数值解, range就是所求数值解的自变量范围。

用Euler方法求初值问题 的数值解。取步长h=0.1,并在一个坐标系中画出数值解与准确解的图形。 解：执行Euler方法程序后，在输入的窗口中按提示分别输入y-2x/y、1、0、1、0.1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果。

y(0.1)=1.1 y(0.2)=1.19182 y(0.3)=1.27744 y(0.4)=1.35821 y(0.5)=1.43513 y(0.6)=1.50897 y(0.7)=1.58034 y(0.8)=1.64978 y(0.9)=1.71778 y(1.)=1.78477

本题的准确解为y(x)= (1 + 2 x)1/2，在一个坐标系中画出数值解与准确解的图形如下。 图中点图是数值解，曲线为准确解。从图中可以知道数值解与准确解比较接近，但还是有误差的。

改进的Euler方法算法 1. 输入函数f(x,y)、初值y0、自变量区间端点a,b步长h 2.计算节点数n和节点x k

改进的Euler方法程序 Clear[x,y,f] f[x_,y_]= Input["函数f(x,y)="] y0=Input["初值y0 ="] a=Input["左端点a="] b=Input["右端点b="] h=Input["步长h="] n=(b-a)/h; For[i=0,i<n,i++, xk=a+i*h; y1=y0+h*f[xk,y0]; y1=y0+h/2*(f[xk,y0]+f[xk+h,y1]); Print["y(",xk+h//N,")=",y1//N]; y0=y1] 注：语句h=Input["步长h="]的步长输入应该用带小数点的数输入，这样可以加快计算速度，否则，由于Mathematica软件本身的原因，可能会出现计算时间过长等计算不出结果的情况。

说明：本程序用改进的Euler公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。 程序中变量说明： f[x,y]: 存放函数f(x,y) y0: 存放初值y0及数值解 a：存放自变量区间左端点 b: 存放自变量区间右端点 n: 存放节点个数 h: 存放节点步长 xk：存放节点xi y1: 存放数值解

Runge-Kutta方法 Runge-Kutta方法算法 1.输入函数f(x,y)、初值y0、自变量区间端点a,b步长h 2.计算节点数n和节点x k 3.用Runge-Kutta r公式求数值解

Runge-Kutta方法程序 Clear[x,y,f] f[x_,y_]= Input["函数f(x,y)="] y0=Input["初值y0 ="] a=Input["左端点a="] b=Input["右端点b="] h=Input["步长h="] n=(b-a)/h; For[i=0,i<10,i++, xk=a+i*h; k1=f[xk,y0]; k2=f[xk+h/2,y0+k1*h/2]; k3=f[xk+h/2,y0+k2*h/2]; k4=f[xk+h,y0+k3*h]; y1=y0+h/6*(k1+2k2+2k3+k4); Print["y(",xk+h//N,")=",y1//N]; y0=y1]

说明：本程序用经典Runge-Kutta公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。 程序中变量说明 f[x,y]: 存放函数f(x,y) y0: 存放初值y0及数值解 a：存放自变量区间左端点 b: 存放自变量区间右端点 n: 存放节点个数 h: 存放节点步长 xk：存放节点xi y1: 存放数值解 注：语句h=Input["步长h="]的步长输入应该用带小数点的数输入，这样可以加快计算速度，否则，由于Mathematica软件本身的原因，可能会出现计算时间过长等计算不出结果的情况。

第六章结束! 谢谢!