数据分析与方程求解PPT课件

1、1,数据分析,2,主要内容,一、多项式计算 二、数据的导入与导出 三、数据的统计分析 四、插值和拟合 五、Mathematica插值 六、Mathematica拟合,3,一、多项式计算,多项式可看做是符号表达式,利用符号运算可进行处理 多项式的向量表示 设有多项式 表示为: (1)多项式的加法 对应系数的加减运算 (2)多项式的乘法 ：conv(P1,P2) (3)多项式除法：Q,r=deconv(P1,P2) 其中Q返回相除的商，r返回余式。Q和r仍是多项式的系数向量 (4)代数多项式求值 ：Y=polyval(P,x) 若x为一数值，则求多项式P在该点的值； 若x为向量或矩阵，则对向量或矩

2、阵中的每个元素求多项式P的值 (5)矩阵多项式求值 ：Y=polyvalm(P,X) X为向方阵，结果为矩阵X的乘与和 (6)多项式的求根 ：x=roots(P) 在复数范围内求多项式P的全部根,4,二、数据的导入与导出,Matlab支持的数据文件 Matlab特制的数据文件：.mat文件 通用数据文件：文本文件、EXCEL文件、数据库文件 等 特制数据文件的打开与保存 保存 右击需保存的变量，通过菜单操作 save命令保存变量到数据文件 打开 通过菜单“File|open” load 命令 通用数据文件 通过“File|Import Data”打开导入向导,5,三、数据的统计分析,求最大值和

3、最小值（max和min） Y,U=max(A) A为向量,将A的最大值存入U,最大值序号存入U A为矩阵,Y表示每列的最大值,U记录每列最大值的行号 U=max(A,B) 由A和B中对应元素的较大者组成矩阵U。 求和、积、均值、中值、累加和、累乘 sum和prod、mean和median、cumsum和cumprod 例：x=1,ones(1,10)*2; y=cumprod(x); sum(y) 求得1+2+22+210 排序 Y,I=sort(A,dim,mode) Y是排序后的矩阵,I记录Y中的元素在A中的位置. dim取1(对列排序);取2(对行排序);默认列排序. mode的取值as

4、cend(升序); descend(降序);默认升序,向量和矩阵,6,四、插值与拟合,引例1:在一天24小时内,从零点开始每间隔2小时测得的环境温度数据分别为(度) 12,9,9,10,18,24,28,27,25,20,18,15,13.请推测中午1点（即13点）时的温度。 引例2:在化学反应中,为研究某化合物的浓度随时间的变化规律,测得一组数据如下： 请给出变化规律,由此推测t为1、1.5、2、2.5、10、10.5、11分钟时的值. 利用已给出的值,计算相关的值 设有一组实验观测数据(xi,yi),i=0,1,n, 如何揭示自变量x与因变量y之间的关系？ 寻找近似的函数关系表达式y=f(

5、x) 常用的方法:插值与拟合,7,1、插值,设测得的n个点的数据为(xi,yi) (i=1,2,n) ,构造一个函数y=g(x),使得在xi(i=1,2,n)有g(xi)=yi。g(x)称为插值函数。 插值函数自变量的个数：一维插值、二维插值、多维插值等 构造插值函数方法:线性插值、多项式插值、样条插值等 一维数据的插值 Y1=interpl(X,Y,X1,mothod) 根据X和Y的值插值，并求插值函数在X1处的值 Y是与X等长的向量,Y1是插值函数在X1处的值 Y是矩阵,X与Y的每一列分别插值，并分别在X1的值 X1是向量或标量,若X1中的元素不在X的范围内,则插值结果为NaN metho

6、d的取值 nearest：最近插值、 linear：线性插值 spline：三次样条插值、cubic：三次插值 二维数据插值 Z1=interp2(X,Y,Z,X1,Y1,method) 根据X、Y和Z的值插值，并求插值函数在X1、Y1处的值 X、Y是向量(X:行,Y:列),Z为函数值(size(Z)=length(Y)length(X) X、Y、Z是矩阵,X与Y同维, Z(i,j)是在点(X(i,j)，Y(i,j)处的值,8,引例1的求解,在一天24小时内,从零点开始每间隔2小时测得的环境温度分别为(度) 12,9,9,10,18,24,28,27,25,20,18,15,13. 请推测中午

7、1点（即13点）时的温度。 解： (1)构造数据 h=0:2:24; T=12,9,9,10,18,24,28,27,25,20,18,15,13; (2)了解数据的分布 绘出数据的图形：plot(h,T,*) (3)选取适当插值方法，计算结果 interp1(h,T,13,spline) interp1(h,T,13,cubic) 插值函数的评价 X=0:24; %选取更多的数据 Y1=interp1(h,T,X,spline); Y2=interp1(h,T,X,cubic); subplot(2,1,1);plot(h,T,:,X,Y1);title(spline) subplot(2,

8、1,2);plot(h,T,:,X,Y2,r);title(cubic,9,例：某实验对一根长为10米的钢轨进行热源的温度传播测试。如下表，其中x表示测量点，h表示测量时间，T表示测得的温度。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度。 解： x=0:2.5:10; h=0,30,60; T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41; x1=0:10; h1=0:20:60; T1=interp2(x,h,T,x1,h1) 等价于 X,H=meshgrid(0:2.5:10,0:30:60); T=95,14,0,0,0;88,48,32

9、,12,6;67,64,54,48,41; X1,H1=meshgrid(0:10,0:20:60); T1=interp2(X,H,T,X1,H1,10,2、拟合,根据一组数据(xi,yi) (i=1,2,n) ，要求确定一个函数y=f(x)，使这些数据点与曲线y=f(x)总体来说尽量接近，称为曲线拟合。 拟合的原理： 最小二乘法 最小 若拟合函数f(x)是一个多项式,称为多项式拟合. Matlab中多项式拟合命令： P=polyfit(X,Y,m) 求数据X与Y的m阶拟合多项式 P为拟合多项式的系数,长为m+1的向量 X与Y是等长的向量 通过polyval求拟合函数在自变量处的值,11,引

10、例2的求解,引例2:在化学反应中,为研究某化合物的浓度随时间的变化规律,测得一组数据如下： 请给出变化规律,推测t为1、1.5、2、2.5、10、10.5、11分钟时的值. 解： (1)构造数据 t=1:10; y=4, 6.4, 8, 8.4,9.28,9.5,9.7, 9.86, 10, 10.2; (2)选取拟合多项式的阶 plot(t,y,o) (3)构造拟合多项式 p=polyfit(t,y,2) p = -0.1109 1.7922 2.9467 (4)计算相关的值 t1=1:0.5:11; y1=polyval(p,t1,12,拟合函数的评价 选取不同的阶做拟合 根据散点图直观判

11、断 计算最小二乘指标 例 t=1:10; y=4, 6.4, 8, 8.4,9.28,9.5,9.7, 9.86, 10, 10.2; p=polyfit(t,y,2); y2=polyval(p,t); plot(t,y,:o,t,y2,-*) sum(y-y2).2) 类似做3次、4次多项式拟合，从中进行比较,13,例：在彩色显影中，由经验得知，形成染料的光学密度与析出银的光学密度由公式 确定，试验测得如下一批数据： 求y关于x的拟合函数。 解：由给定的经验公式来求拟合函数，无法直接用polyfit函数，我们换一种思路，将经验公式两边取对数，得 令Y=lny，X=1/x 只需求Y关于X的线

12、性拟合 x=0.05,0.06,0.07,0.10,0.14,0.20,0.25,0.31,0.38,0.43,0.47; y=0.10,0.14,0.23,0.37,0.59,0.79,1.00,1.12,1.19,1.25,1.29; X=1./x; Y=log(y) P=polyfit(X,Y,1) P = -0.1459 0.5476 则Y=-0.1459X+0.5476 y=exp(Y)=exp(0.5476)e(-0.1459X)=exp(0.5476)e(-0.1459/x),整理即可,14,五、Mathematica插值,构造插值对象的函数： Interpolationdata

13、,InterpolationOrdern 功能： 对数据data进行插值,并可设置插值多项式的次数n, 默认值为3。 注: 生成一个InterpolatingFunction插值范围，目标 所得目标为近似函数,不显示所构造的函数。 数据表示：x0,f0,x1,f1,xn,fn (平面点的坐标) 引例1求解解： data=0,12,2,9,4,9,6,10,8,18,10,24,12,28,14,27,16,25,18,20,20,18,22,15,24,13 ListPlotdata f=Interpolationdata f(13,15,六、Mathematica拟合,格式： Fit数据，拟

14、合函数的基，变量 功能： 用数据data，按给定的变量和拟合函数的基构造拟合函数。 常用的几种格式： Fitdata,1,x,x 用数据data作线性拟合函数a+bx Fitdata,Tablexi,i,0,n,x 作n次多项式拟合 引例2的求解： 构建数据表： data=1,4,2,6.4,3,8.0,4,8.4,5,9.28,6,9.5,7,9.7,8,9.86,9,10,10,10.2 选取拟合的基。 画图，分别作3次、4次、5次、6次多项式拟合，求出其误差，比较、判断几次拟合较为合适,16,插值与拟合 插值要求函数在每一个观测点处一定要满足yi=f(xi). 拟合主要考虑到观测数据受随

15、机误差的影响,寻求整体误差最小、较好反映观测数据的近似函数. 插值与拟合一般过程： 构造数据 确定插值多项式的阶（可采用经验公式或画图） 选取命令插值 求出并分析所得结果 数据拟合的一般过程： 考察数据来源,确定有无经验公式； 若有经验公式,则可直接进行拟合（必要时可对数据作些变换处理）; 若无经验公式,选择拟合函数的类型（画图） 求出拟合函数 从数据、图形等途径对拟合结果分析; 确定结果或调整拟合函数,17,方程求解,18,主要内容,符号方程的求解 符号代数方程求解 符号常微分方程求解 线性方程组求解 利用左除运算符 利用矩阵的分解 非线性方程求解 单个自变量的非线性方程求解 非线性方程组求

16、解 微分方程求解 Mathematica方程的求解,19,一、符号方程的求解,符号代数方程求解 代数方程是指未涉及微积分运算的方程 求解命令 solve(s) 求解符号代数方程，求解变量是默认变量。 solve(s1,s2,.,sn,v1,v2,.,vn) 求解符号表达式s1,s2,.,sn所组成的代数方程组的解，求解变量分别为v1,v2,.,vn。 例： solve (x2-a*x-b=0) solve (x2-a*x-b) solve (x2-a*x-b=0 , x ) syms x a b f=x2-a*x-b; solve(f,x); solve(f) syms x y z x y z

17、=solve(x+2*y-z=27,x+z=3,x2+3*y2=12,x,y,z) f=x5-a*x-b=0;solve(f,x) %无法求出解析解 f=x5-x-3;solve(f) %求近似解,20,符号常微分方程求解 方程的表示： Dy表示y，D2y表示y 求解命令 dsolve(e,c,v) 求解微分方程e在初始条件c下的特解,参数v描述方程中的自变量. 例： (1)求 当y(0)=1时的特解。 y=dsolve(Dy=2*y*y2,y(0)=1,x) (2)求 的通解。 x,y=dsolve(Dx=4*x-2*y,Dy=2*x-y,t,21,二、线性方程组求解,若有线性方程组 利用左

18、除运算符 A=2,1,-1,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; x=Ab 利用矩阵的分解 LU分解：将矩阵表示为一个下三角矩阵与一个上三角矩阵的乘积。 L,U=lu(X): 产生L和U ,使得X=LU。 A=2,1,-1,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; L,U=lu(A); x=U(Lb) QR分解:是将矩阵分解为一个正交矩阵和一个上三角矩阵的乘积. ： Q,R=qr(X): 产生Q和R,使得X=QR,22,三、非线性方程求解(近似解,例：考察函数 作图观察其根的分布 尝试用solve求

19、解 单个自变量的非线性方程求解 z=fzero(fun, x0) z=fzero(fname,x0) 在x0附近寻找函数fun的近似根(fname是待求根的函数名) 求上述函数在区间-3,1内的根 （1）画图确定根的范围 f=abs(x*sin(x)-exp(x)-1; fplot(f,-3,1) 观察图形得出在x=0、-1.2、-2.6附近有根 （2）调用fzero函数分别求出其根 fzero(f,0) fzero(f,-1.2) fzero(f,-2.6,23,上述问题也可这样求解 (1)定义函数文件 funx.m, function y=funx(x) y=abs(x.*sin(x)-e

20、xp(x)-1; (2)画图确定根的位置 fplot(funx,-3,1) %funx为函数funx的句炳 观察图形得到出在x=0、-1.2、-2.6附近有根 (3)调用fzero函数分别求出其根 fzero(funx,0) fzero(funx,-1.2) fzero(funx,-2.6) 上述问题还可这样求解 (1)定义内联（匿名）函数 ff=(x)abs(x.*sin(x)-exp(x)-1; (2)画图确定根的位置 fplot(ff,-3,1) (3)调用fzero函数分别求出其根 fzero(ff,-2.5) fzero(ff,-1.5) fzero(ff,0,24,非线性方程组求解

21、,求解方程组 初始值为-5,5 MATLAB的优化工具箱（Optimization Toolbox） 对于非线性方程组F(x)=0,求解命令为： X = fsolve(fname, X0, options) 解 (1)定义函数文件myfun.m function f = myFun(x) f(1) = 2*x(1) - x(2) - exp(-x(1); f(2)= -x(1) + 2*x(2) - exp(-x(2); (2) 提交给求解方程的函数 -5,-5 %初始值 op=optimset(display,off); %修改参数选项,中间结果不显示 x=fsolve(myfun, x0,

22、 op) x = 0.5671 0.5671,25,四、微分方程求解,求解微分方程： MATLAB提供了多个求常微分方程数值解的函数,格式为: t,y=solver(fname,tspan,y0) solver为求常微分方程数值解的函数 t和y分别给出时间向量和相应的状态向量。 fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。 tspan为求解区间形式为t0,tf,y0是初始状态列向量。 解： 建立函数文件funt.m function yp=funt(t,y) yp=(y2-t-2)/4/(t+1); 求解微分方程 t,y=ode23(funt,0,10,2); %求

23、数值解,26,五、Mathematica方程的求解,方程的表示 方程中的等号（=）应用逻辑等号（= =） 根的表示：逻辑量、转换规则 多项式方程求解 Solve方程或方程组,变量列表 NSolve方程或方程组,变量列表 注： Solve主要是处理多项式方程，尽可能给出精确解，若给不出显示公式解，给出隐式公式解。 例：求方程x5-x+11=0的根。 NSolve求出的是近似数值解(可求任意多项式方程) 方程组的表示： 方式1：lhs1= =rhs1, lhs2= =rhs2 方式2：lhs1= =rhs1 g=diff(f) %求导函数 sin(x)+x*cos(x)-exp(x) fzero(

24、sin(x)+x*cos(x)-exp(x),-2) 结果： -2.0745 故在x-2.0745处函数取得极大值。 上述问题也求解如下： syms x f=x*sin(x)-exp(x); fplot(char(f),-3,2); %将符号表达式转换为字符串 x0=fzero(char(diff(f),-2); y0=subs(f,x,x0); %将符号表达式f中的x替换为x0,求得极大值. x0,y0 -2.0745,1.6912,31,Matlab专门提供了求极小值的函数，格式为： (1) x,fval=fminbnd(fname,x1,x2,option) 求一元函数在区间(x1,x2

25、)中的极小值点x和极小值。 (2) x,fval=fminsearch(fname,x0,option) 用单纯形法求多元函数在x0附近的极小值点x和极小值. (3) x,fval=fminunc(fname,x0,option) 用拟牛顿法求多元函数在x0附近的极小值点x和极小值. 例:求y=x*sin(x)-ex在区间-3,2内的极值点。 (1)建立函数文件myfin.m,命令如下： function y=myfin(x) y=-(x.*sin(x)-exp(x); (2)调用fminbnd函数求极大值,命令如下： x,fval=fminbnd(myfin,-3,-1.5) 结果：x = -2.0745，fval = -1.6912 y=x*sin(x)-ex在x = -2.0745处取得极大值1.6912。 等价于 f=(x)-(x*sin(x)-exp(x); x,fval=fminbnd(f,-3,-1.5,32,二、有约束问题求解,模型： 求一组x ( x=x1,x2,xnT)使得目标函数f(x)为最小,且x满足约束条件G(x)0. 约束条件可表示为： 线性不等式约束：A