数学建模MATLAB教案

1、数学建模 MATLAB 教案1 / 61三种插值方法 拉格朗日多项式插值 构造基函数l(X)(x Xo)L (x XiJ(x XiJL (x Xn)i(X Xo)L (XiXi i)(XiXiJL (XiXn)插值多项式nLn(x)yi(x)i 0分段线性插值将每两个相邻的节点用直线连起来，即在每个小区间上是线性函数。有现成命令。 三次样条插值一根有弹性的细长木条固定在节点上，其他地方自然弯曲，如此称为样条曲线。 普遍使用的样条函数是分段三次多项式：在每个小区间上是三次多项式，在大区间上二阶导数连续，通过全部节点。有现成命令。例子1对 y 亠，5x5，用 11 个等分节点作上述三种插值，用 2

2、1 个等分插值点1 x作图。比较结果，spli ne 插值最好。2.数据拟合2.1 多项式拟合指令方法x=1 2 3 4 5 6 7 8 9;y=9 7 6 3 -1 2 5 7 20;P=polyfit(x,y,3);xi=0:.2:10;yi=polyval(P,xi);plot(xi,yi, x,y,r*);图形窗口方法x=1 2 3 4 5 6 7 8 9; y=9 7 6 3 -1 2 5 7 20; plot(x,y,r*);2.2 指定函数类型拟合例子某次阻尼振荡实验中测得数据点x=0;0.4;1.2;2;2.8;3.6;4.4;5.2;6;7.2;8;9.2;10.4;11.6

3、;12.4;13.6;14.4;15;y=1;0.85;0.29;-0.27;-0.53;-0.4;-0.12;0.17;0.28;0.15;-0.03;-0.15;-0.071;0.059;0.08;0.032;-0.015;-0.02;f=fittype(a*cos(k*t)*exp(w*t),i ndepe nden t,t,coefficie nts,a,k,w);cfun=fit(x,y,f)xi=0:.1:20;yi=cfu n( xi);plot(x,y,r*,xi,yi,b-);补充：三维图形例子数学建模 MATLAB 教案2 / 6作曲面z2 2f (x, y)的图形，zsi

4、y,x2, 7.5 x 7.5, 7.5 y 7.5Jx2y2作螺旋线x sint,y cost,z t练习矩阵运算：+加法；-减法；转置；*乘法；八乘幕；左除；/右除。.*乘法；A乘幕；左除；./右除。矩阵函数：zeros；ones eye；rand; randn。图形：线型：-实线；：点线；-.虚点线；-波折线；圆点；+加号；*星号；x x 形; 小圆。颜色：y 黄；r 红；g 绿；b 蓝；w 白；k 黑；m 紫；c 青。例子车灯光源投影区域的绘制(CUMCM2002A)3.数值积分矩形公式Lnn 1hfk,hk 0fk,h梯形公式(相当于将两矩阵公式平均，也相当于用分段线性插值函数作为近

5、似)n 1hTnhfk-(fo+fn),hk 12辛普森公式(抛物线法)m 1m 1( f0+f2m4f2k 12 f2k), h3k 0k 1b a2m例子02si nxdx例子卫星轨道长度人造地球卫星轨道可视为平面上的椭圆。我国第一颗人造地球卫星近地点距地球 表面 439km,远地点距地球表面 2384km,地球半径为 6371km,求该卫星的轨道 长度。4.数值微分f(a)f(a)f(a)f(af(a h)h) f f前差公式 hf ff(af(a h)h)后差公式 hf(af(a h)h) f(af(a h)h)中心差商2h数学建模 MATLAB 教案3 / 6数学建模 MATLAB

6、教案4 / 6前差形式数值微分，现成命令：diff(x)。输入 x 是 n 维数组，输出为 n-1 维数组5常微分方程数值解设yf (X, y), y(Xo)乂，其中 f 适当光滑，对 y 满足 Lipschitz 条件，以保证解存在且唯一。在一系列离散点x0 x-ix2L xnL上求y(xn)的近似值yn(n 1,2,L )，通常向前欧拉公式在小区间Xn,Xm上用差商y(Xy(Xnl) ) y(Xy(Xn) )代替方程左端的导数，方程右端f(X,y)h中的X取小区间召，召i的左端点Xn，可得y(Xn i) y(Xn) hf(Xn,y(Xn)从Xo出发，由初值y(Xo) yo，得到y(xj的近

7、似值y-为(以下用 代替)yiyohf (Xo, yo)再以y-作为y(xi)的近似值，得到y(X2)的近似值y为y2yihf (Xi, yi)继续下去，一般地yn iynhf (Xn,yn),n0,i,L龙格库塔方法按照微分中值定理有yxyx ww y(xnh),o ih注意到方程y f (x, y)就有y(Xn i) y(Xn) hf (Xnh, y(Xnh)记K f (Xnh, y(Xnh)，称为区间Xn,Xni上的平均斜率。向前欧拉公式简单地取f (xn, yn)为 K，精度很低。龙格库塔方法的基本思想：在区间Xn,Xni内多取几个点，将它们的斜率加权平 均作为 K 取 2 个点，2

8、阶龙格库塔公式。取 4 个点，4 阶龙格库塔公式。有现成的命令：X2Xi,X3X2,L ,xnXn 1。取等步长 h，即Xnx0nh(n 1,2,L )。数学建模 MATLAB 教案5 / 6t,x=ode23(f,ts,x0,optio ns)t,x=ode45(f,ts,x0,optio ns)其中，f 是由待解方程写成的 m 文件名；ts 为自变量的取值；x0 为函数的初值； options 用于设定误差限(可以缺省，缺省时设定为相对误差103，绝对误差106)， 设定命令：options=odeset(reltol,rt,abstol,at) rt， at 分别为设定的相对误差和绝 对

9、误差。t，x 为输出的自变量和函数值。例子单摆运动方程ml mg sin，初始条件(0)0,牌0) 0在0不大的条件下，可将方程中的 sin 近似为，于是得到线性常系数微分方程 皱g0，容易算出其在初始条件下的解为(t)0cos t,、g。当较大时，若仍用 近似 sin ，误差太大了。试用数值方法在等于10和30o两种情况下求解(设 l 25cm)，画出(t)的图形，并与近似解的结果比较。先将它化为方程组。令Xi,X2&，则方程化为& x2,x2gsin，初始条件为x1(0)为0,x2(0)0，其中g 9.8， l 25， 细 为10o0.1745弧度和30o0.5236弧度两

10、种情况。对于近似解，周期 T 2 ,学 2 ,10(s)可以看出，初始角度为10o时精确做值)解与近似解相差不大，而初始角度为30o时，随着时间的增加二者差别就很大了。食饵-捕食者模型 食饵甲和捕食者乙在时刻t的数量分别记作x(t)和y(t)当甲独立生存时它的(相对)增长率为r，即& rx，而乙的存在使甲的增长率减小， 设减小的程度与乙的数量成正比，于是x(t)满足2 x(r ay) rx axy，比例系数a反映捕食者掠取食饵的能力。捕食者离开食饵无法生存，设乙独自存在时死亡率为d，即& dy，而甲的存在使乙的死亡率降低，设这个作用与甲的数量成正比，于是y(t)满足&

11、y( d bx) dy bxy，比例系数 b 反映食饵对捕食者的供养能力。设初始数量分别为x(0) Xo, y(0)yo。数学建模 MATLAB 教案6 / 6设r 1,d0.5,a 0.1,b 0.02, Xo25, y。2，求数值解，画出函数x(t)和y(t)的图形以及相图(x,y)6.线性规划LINPROGLin ear program ming.X=LINPROG(f,A,b) solves the linear programming problem:min f*xxsubject to: A*x = b例子min z2为3x2x34x22x38st3x12x26人丛， 0c=2;3

12、;1;a=1,4,2；3,2,0;b=8;6;x,y=li nprog(c,-a,-b,zeros(3,1)c=2;3;-5;a=-2,5,-1;b=-10;aeq=1,1,1;beq=7;x=li nprog(-c,a,b,aeq,beq,zeros(3,1)value=c*x7.优化工具箱help optim 查看优化工具箱的帮助信息基本程序名：fmi n, fmi nu,fmi ns,leastsq,curvefit,co nstr,lp,qp,mi ni max, nnl s,co nls,fzero,fsolve等。例子2 2求解min?a(本题精确结果为 x=y=0)用下面一组数据拟合c(t) rekt中的系数 r 和 k。(也可以用指定函数类型拟合)max z2%3x25x3% x?X37st 25x2X310 xn%2压0数学建模 MATLAB 教案7 / 6t=0.25 0.5 1 1.5 2 3 4 6 8;c=19.21 18.15 15.36 14