5数值计算.docx
- 文档编号:3542301
- 上传时间:2022-11-23
- 格式:DOCX
- 页数:13
- 大小:50.37KB
5数值计算.docx
《5数值计算.docx》由会员分享,可在线阅读,更多相关《5数值计算.docx(13页珍藏版)》请在冰豆网上搜索。
5数值计算
第六章数值计算命令与例题
数值问题由一组已知数据来求出一组结果数据,使得这两组数据之间满足预先制定的某种关系的问题,数值计算就是通过对数学问题解有效的近似处理去获得比较满意解的计算方法。
这里只介绍与数值计算中最基本方法的一些命令与例题。
6.1求近似函数
在生产和实验中,人们经常遇到需要通过某个未知(或复杂)的函数f(x)在有限个给定点的函数值:
{xi,yi},i=1,2,….,n,这里f(xi)=yi去获得函数f(x)的近似函数(x),以便于处理涉及到f(x)的计算或其他处理问题,这种求近似函数(x)的方法主要有拟合方法和插值方法。
一般,曲线拟合是常用的拟合方法,而插值方法常用的有多项式插值与分段插值方法,在Mathematica中它们都有相应的命令来获得相应的近似函数。
6.1.1曲线拟合
曲线拟合主要用来求一元近似函数,它是根据一组函数值(通常是测量值或试验数据)去寻找描述相应函数关系的一种常用数据逼近方法,是根据最小二乘原理的意义下获得近似函数的,此近似函数具有在数据点处的误差平方和最小的特点。
通过拟合得到的函数关系通常称为经验公式,该函数不必经过任何的数据点,是实际函数的一种近似函数。
曲线拟合结果的好坏与选用什么类型的函数作为拟和曲线类型有关。
记函数集合:
M=Span[0,1,2,…,m]={(x)|(x)=a00(x)+a11(x)+…+amm(x),aiR}
称集合M为函数0,1,2,…,m张成的空间,m+1个函数0(x),1(x),2(x),…,m(x)称为拟合基函数集合,它们都是已知的函数。
拟合基函数集合不同对应不同的拟合函数类,如果要用曲线拟合方法得到经验公式,通常是先画出所给数据的散点图,再根据散点图选择合适的拟和函数类型做曲线拟合。
应该注意的是:
选择不同的拟合函数类会得到不同的拟合函数,拟合函数的好坏依赖于所选的拟合函数类,如果你知道的曲线类型多,就容易获得满意的拟合函数。
利用Mathematica提供的曲线拟合命令可以方便快速地尝试各种不同的拟合函数类型可以帮你找到满意的近似函数。
Mathematica曲线拟合的一般形式为:
Fit[{数据点集合},{拟合基函数集合},自变量名]
具体的拟合命令有:
命令形式1:
Fit[{{x1,y1},{x2,y2},...,{xn,yn}},{0,1,2,…,m},x]
功能:
根据数据点集{{x1,y1},{x2,y2},...,{xn,yn}}求出具有拟合函数为
(x)=a00(x)+a11(x)+…+amm(x)
形式的近似函数(x)
命令形式2:
Fit[{y1,y2,...,yn},{0,1,2,…,m},x]
功能:
根据数据点集{{1,y1},{2,y2},...,{n,yn}}求出具有拟合函数为
(x)=a00(x)+a11(x)+…+amm(x)
形式的近似函数(x)
命令形式3:
Fit[{{x1,y1},{x2,y2},...,{xn,yn}},Table[xi,{i,0,m}],x]
功能:
根据数据点集{{x1,y1},{x2,y2},...,{xn,yn}}求出拟合函数为m次多项式的近似函数(x)=a0+a1x+a2x2+…+amxm
注意:
命令1和命令2中的0,1,2,…,m可以是任何具体的已知函数。
例题
例1.已知数表{{-3,-1.2},{-1,1.3},{0,1.5},{1,1.9},{3,2}},求二次拟合曲线方程。
解:
二次拟合曲线就是拟合函数为二次多项式类型,它的拟合基函数为
{1,x,x2},因此Mathematica命令为:
In[1]:
=Fit[{{-3,-1.2},{-1,1.3},{0,1.5},{1,1.9},{3,2}},{1,x,x2},x]
Out[1]=1.65238+0.51x-0.138095x2
例2.产生函数值在0到5.0之内波动的20个随机数值,分别用二次,三次和七次拟合曲线来拟合这组数据,说明这三个拟合函数哪个效果更好?
解:
Mathematica命令:
In[2]:
=t=Table[5*Random[],{20}];*产生函数值在0到5.0之内的20个随机数值
Out[2]={1.69283,2.02028,0.446517,3.00438,3.50901,3.53736,4.84185,2.37896,3.96983,3.40952,1.93415,1.42763,2.48143,2.03451,0.708301,2.60303,4.08826,0.203527,3.36201,0.08036}
In[3]:
=g=ListPlot[t,PlotStyle->PointSize[0.01]]*将散点图图形文件存放在变量g中
In[4]:
=f1=Fit[t,{1,x,x^2},x]*将二次拟合函数存放在变量f1中
Out[4]=1.53945+0.31661x-0.0172625x2*获得的二次拟合函数
In[5]:
=s1=Plot[f1,{x,1,20}]*将二次拟合函数图形文件存放在变量s1中
In[6]:
=Show[s1,g]*将散点图和二次拟合函数图画在同一个坐标系中
In[7]:
=f2=Fit[t,{1,x,x^2,x^3},x]*将三次拟合函数存放在变量f2中
Out[7]=0.430077+0.882677x-0.0830357x2+0.00208804x3*获得的三次拟合函数
In[8]:
=s2=Plot[f2,{x,1,20}]*将三次拟合函数图形文件存放在变量s2中
In[9]:
=Show[s2,g]*将散点图和三次拟合函数图画在同一个坐标系中
In[10]:
=f3=Fit[t,Table[x^i,{i,0,7}],x]*将七次拟合函数存放在变量f3中
Out[10]=5.09229-5.01425x+2.15792x2-0.375783x3+0.0332493x4
-0.00163526x5+0.0000442756x6-5.3285410-7x7
*获得的七次拟合函数
In[11]:
=s3=Plot[f3,{x,1,20}]*将七次拟合函数图形文件存放在变量s3中
In[12]:
=Show[s3,g]*将散点图和七次拟合函数图画在同一个坐标系中
通过比较,可以知道,显然七次拟合函数效果较好
例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},试用合适的函数进行拟合
解:
Mathematica命令:
In[13]:
=t={2.89229,2.86323,0.473147,-2.25209,-2.87003,
-0.835768,1.97187,2.96841,1.23648,-1.63202};
In[14]:
=g=ListPlot[t,PlotStyle->PointSize[0.04]]*画出散点图
In[15]:
=f=Fit[t,{1,Sin[x]},x]*由散点图选(x)=a0+a1sin(x)作为拟合函数类做尝试
Out[15]=0.0482784+3.07027Sin[x]*获得拟合函数
In[16]:
=s=Plot[f,{x,1,10}]*将拟合函数图形文件存放在变量s中
In[17]:
=Show[s,g]*将散点图和拟合函数图画在同一个坐标系中
Out[17]=
-Graphics-
从图中可以知道,拟合函数效果很好。
6.1.2函数插值
插值法由实验或测量的方法得到所求函数y=f(x)在互异点x0,x1,...,xn处的值y0,y1,…,yn,构造一个简单函数(x)作为函数y=f(x)的近似表达式
y=f(x)(x)
使(x0)=y0,(x1)=y1,,(xn)=yn,(x)称为插值函数,它常取为多项式或分段多项式。
与曲线拟合函数不同的是插值函数(x)满足条件(x0)=y0,(x1)=y1,,(xn)=yn。
●多项式插值
多项式插值是在给定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)
例题
例3.已知数表{{2,1.41421},{2.1,1.44914},{2.2,1.48324},{2.4,1.54919}},求其插值多项式。
解:
因为是4个点,故它的插值多项式次数是3次
因此Mathematica命令为:
In[18]:
=InterpolatingPolynomial[{{2,1.41421},{2.1,1.44914},{2.2,1.48324},{2.4,1.54919}},x]
Out[18]=1.41421+(0.3493+(-0.0415+0.01(-2.2+x))(-2.1+x))(-2+x)*插值多项式次
In[19]:
=Simplify[%]
Out[19]=0.44891+0.65165x-0.1045x2+0.01x3*化简后的插值多项式次
例4.观察函数f(x)=(1+x2)-1在区间[-4,4]内选取9个等距插值节点所得8次插值函数与函数f(x)=(1+x2)-1关系,说明高次插值函数效果不好
解:
Mathematica命令:
In[20]:
=u=Table[{x,(1+x^2)^-1},{x,-4,4}];*产生函数f(x)=(1+x2)-1的9个等距插值节点In[21]:
=g=ListPlot[u,PlotStyle->PointSize[0.04]]*将散点图图形文件存放在变量g中
In[22]:
=s=InterpolatingPolynomial[u,x];*将插值函数存放在变量s中
In[23]:
=(1+x^2)^-1/.x->3.7*计算被插值函数在x=3.7的函数值
Out[23]=0.0680735
In[24]:
=s/.x->3.7*用插值函数计算在x=3.7的函数值
Out[24]=-0.662229*可以看到在x=3.7误差很大
In[25]:
=t=Plot[{s,(1+x^2)^-1},{x,-4,4},
PlotStyle->{{Thickness[0.005]},{Thickness[0.006]}}]
*将插值函数s与f(x)画在一起的图形文件存放在变量t中
In[26]:
=Show[t,g]*将散点图,插值函数s,f(x)画在一个坐标系中
Out[26]:
=-Graphics-
图中细的曲线是插值函数,粗的曲线是原来函数(被插函数),从中可以看到在靠近两端误差很大,故高次插值函数误差较大,这也是为什么实用中一般不使用多项式插值的原因。
●分段多项式插值
分段多项式插值是在给定n个数据点:
{xi,yi},i=1,2,….,n,a=x0 (x)在区间[a,b]上连续;(xi)=yi (x)在每个子区间[xi,xi+1](i=0,1,2,,n-1)上是次数不超过m的多项式 这样的(x)称为是f(x)在[a,b]上的分段m次插值多项式函数。 常用的是m=3的分段3次插值多项式。 分段多项式插值不会出现高次多项式插值引起的误差增大的问题,实际和理论可以证明,随着插值点的增多,与多项式插值不同的是,分段插值函数的误差会越来越小。 因此,实用中常使用分段多项式插值作为插值函数。 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) 注意: 用Mathematica命令求出的分段插值多项式没有给出具体的分段函数表达式,而是用 “InterpolatingFunction[{x1,xn},<>]”作为所求的分段插值函数,通常可以用 变量=Interpolation[{数据点集合}] 把所得的分段多项式插值存放在变量中,如果要计算插值函数在某一点的如p点的值,只要输入“变量[p]”即可,如果要表示这个插值函数应该用“变量[x]”。 例题 例5.用分段多项式插值重做例4的插值问题,并验证,随着插值点的增多,分段插值函数的误差会越来越小。 解: Mathematica命令: In[27]: =r=Interpolation[Table[{x,(1+x^2)^-1},{x,-4,4}]]*将插值函数存放在变量r中 Out[27]=InterpolatingFunction[{-4,4},<>]*形成在[-4,4]上的分段插值函数 In[28]: =(1+x^2)^-1/.x->3.7*计算被插值函数在x=3.7的函数值 Out[28]=0.0680735 In[29]: =r[3.7]*用分段插值函数计算在x=3.7的函数值Out[29]=0.0734*可以看到在x=3.7误差较小 In[30]: =d=Table[{x,(1+x^2)^-1},{x,-4,4,0.5}];*加密插值节点并将数据点存在变量d中 In[31]: =r1=Interpolation[d]*将插值函数存放在变量r1中 Out[31]=InterpolatingFunction[{-4,4.},<>]*形成在[-4,4]上的分段插值函数 In[32]: =r1[3.7]*用插值函数计算在x=3.7的函数值 Out[32]=0.0681761*可以看到在x=3.7误差更小 In[33]: =Plot[{r1[x],(1+x^2)^-1},{x,-4,4}]*将插值函数r(x)与f(x)画在同一个坐标显示 Out[33]=-Graphics- 图中可以看出分段插值函数与被插函数误差在整个插值区间[-4,4]都很小,两个图形几乎相同,插值效果非常好。 结果显示,随着插值点的增多,分段插值函数的误差越来越小。 6.2解常微分方程 由自变量、未知函数以及未知函数的导数(或微分)之间的一个(或一组)关系式称为微分方程(或微分方程组)。 未知函数都是同一个自变量的微分方程(或微分方程组)称常微分方程(或常微分方程组)。 微分方程(或微分方程组)中出现的未知函数最高阶导数的阶数,称为微分方程(或微分方程组)的阶数。 求出常微分方程(或微分方程组)的未知函数(或未知函数组),称为解常微分方程。 含有任意常数的解称为通解,否则称为特解。 n阶常微分方程的一般形式为: 利用Mathematica可以象高等数学一样求出常微分方程的解析解(准确解),如果求不出解析解Mathematica可以求出微分方程的数值解(即近似解),有关微分方程的数值解的知识可以参看文献[3]的常微分方程初值问题的数值解法内容。 Mathematica解常微分方程的命令有求常微分方程(组)的解析解命令和求常微分方程(组)的数值解命令。 6.2.1求常微分方程(组)的解析解 命令形式1: DSolve[eqn,y[x],x] 功能: 求出常微分方程eqn的未知函数y[x]的解析通解。 命令形式2: DSolve[{eqn1,eqn2,...},{y1[x],y2[x],...},x] 功能: 求出常微分方程组{eqn1,eqn2,...}的所有未知函数{y1[x],y2[x],...}的解析通解。 注意: 1.每个常微分方程的中的等号输入时应该用两个等号代替 2.未知函数及未知函数的导数要用“[自变量名]”指示未知函数的自变量 3.常微分方程组和未知函数组应该用大括号括起来 4.命令形式2可以用来求常微分方程的特解 5.eqn表示常微分方程,x是自变量名,y[x]表示未知函数,自变量名和未知函数可以是其他的变量名 例题 例6.求常微分方程y=2ax的通解,a为常数 解: Mathematica命令: In[1]: =DSolve[y'[x]==2ax,y[x],x] Out[1]={{y[x]->ax2+C[1]}} 即得本问题的通解 例7.求常微分方程y+y=0的通解 解: Mathematica命令: In[2]: =DSolve[y''[x]+y[x]==0,y[x],x] Out[2]={{y[x]->C[2]Cos[x]-C[1]Sin[x]}} 即得本问题的通解: C1,C2是任意常数 例8. 求常微分方程 的特解。 解: Mathematica命令: In[3]: =DSolve[{y'[x]==x/y[x]+y[x]/x,y[1]==2},y[x],x] Out[3]={{y[x]->Sqrt[x2(4+2Log[x])]}} 本问题的特解为: 下面的操作可以检验所求特解的正确性: In[4]: =y[x_]: =Sqrt[x^2*(4+2Log[x])]*将所求特解定义为一个函数 In[5]: =y[1]*检验初始条件y[1]==2 Out[5]=2*初始条件成立 In[6]: =D[y[x],x]-x/y[x]-y[x]/x;*将微分方程左端项与右端项相减 In[7]: =Simplify[%]*化简计算结果 Out[7]=0*说明所求特解正确 例9.求常微分方程组y(t)=x(t),x(t)=y(t)的通解 解: Mathematica命令: In[8]: =DSolve[{y'[t]==x[t],x'[t]==y[t]},{x[t],y[t]},t] C[1]+E2tC[1]-C[2]+E2tC[2] Out[8]={{x[t]->----------------------------------------, 2Et -C[1]+E2tC[1]+C[2]+E2tC[2] y[t]->----------------------------------------}} 2Et 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,...的数值解。 注意: 1)求常微分方程(组)的数值解命令中的常微分方程(组)中不能出现除未知函数及其导数和自变量之外的任何其他未知参数字母。 2)由NDSolve命令得到的解是以{{未知函数名->InterpolatingFunction[range,<>]}}的形式给出的,其中的InterpolatingFunction[range,<>]是所求的插值函数表示的数值解,range就是所求数值解的自变量范围。 3)InterpolatingFunction[range,<>]是不带自变量的函数名,其函数表示为InterpolatingFunction[range,<>][自变量名]。 4)eqns表示常微分方程(组)及初值条件,x是自变量名,y或y1,y2,…表示未知函数名,注意未知函数名是指不必用方括号指示其自变量的函数名,如函数f(x)的函数名为f,正弦函数sinx的函数名为sin。 自变量名和未知函数可以是其他的变量名。 例题 例10.求微分方程初值问题y=x+y,y(0)=1在区间[0,0.5]内的数值解 解: Mathematica命令: In[1]: =NDSolve[{y'[x]==x+y[x],y[0]==1},y,{x,0,0.5}] Out[1]={{y->InterpolatingFunction[{0.,0.5},<>]}} 上面显示的是所求数值解的替换形式,为得到本问题数值解,再键入: In[2]: =f=y/.%[[1]]*把Out[1]中的InterpolatingFunction[{0.,0.5},<>]赋值给变量f Out[2]=InterpolatingFunction[{0.,0.5},<>] In[3]: =Plot[f[x],{x,0,0.5}]*画出所求未知函数数值解的图形 Out[3]=-Graphics- In[4]: =Table[f[x],{x,0,0.5,0.1}]*计算数值解在0,0.1,0.2,0.3,0.4,0.5的值 Out[4]
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算
![提示](https://static.bdocx.com/images/bang_tan.gif)