5数值计算Word下载.docx
《5数值计算Word下载.docx》由会员分享，可在线阅读，更多相关《5数值计算Word下载.docx（13页珍藏版）》请在冰豆网上搜索。
Fit[{y1,y2,...,yn},{0,1,2,…,m},x]
根据数据点集{{1,y1},{2,y2},...,{n,yn}}求出具有拟合函数为
命令形式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},试用合适的函数进行拟合
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）
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关系,说明高次插值函数效果不好
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<
x1<
x2<
<
xn=b后,求出一个近似函数（x）作为函数y=f（x）的近似函数,（x）满足条件
（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[{数据点集合}]
具体的分段多项式插值命令有:
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]:
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.},<
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求常微分方程（组）的解析解
DSolve[eqn,y[x],x]
求出常微分方程eqn的未知函数y[x]的解析通解。
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为常数
In[1]:
=DSolve[y'
[x]==2ax,y[x],x]
Out[1]={{y[x]->
ax2+C[1]}}
即得本问题的通解
例7.求常微分方程y+y=0的通解
'
[x]+y[x]==0,y[x],x]
Out[2]={{y[x]->
C[2]Cos[x]-C[1]Sin[x]}}
即得本问题的通解:
C1,C2是任意常数
例8.
求常微分方程
的特解。
=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*初始条件成立
=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）的通解
[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]->
----------------------------------------}}
6.2.2求常微分方程（组）的数值解命令
NDSolve[eqns,y,{x,xmin,xmax}]
求出自变量范围为[xmin,xmax]且满足给定常微分方程及初值条件eqns的未知函数y的数值解
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]内的数值解
=NDSolve[{y'
[x]==x+y[x],y[0]==1},y,{x,0,0.5}]
Out[1]={{y->
InterpolatingFunction[{0.,0.5},<
]}}
上面显示的是所求数值解的替换形式,为得到本问题数值解,再键入:
=f=y/.%[[1]]*把Out[1]中的InterpolatingFunction[{0.,0.5},<
]赋值给变量f
Out[2]=InterpolatingFunction[{0.,0.5},<
]
=Plot[f[x],{x,0,0.5}]*画出所求未知函数数值解的图形
Out[3]=-Graphics-
=Table[f[x],{x,0,0.5,0.1}]*计算数值解在0,0.1,0.2,0.3,0.4,0.5的值
Out[4]
升级会员 