建模matlab之三符号计算及其应用.ppt

《建模matlab之三符号计算及其应用.ppt》由会员分享，可在线阅读，更多相关《建模matlab之三符号计算及其应用.ppt（122页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、MMATLABATLAB 软件及其应用王林君江苏大学理学院 Application of Matlab Language1MATLAB符号运算符号运算（Symbolic）2Matlab符号运算介绍q Matlab 符号运算是通过符号运算是通过符号数学工具箱符号数学工具箱（Symbolic Math Toolbox）来实现的。）来实现的。Matlab 符号数学工具箱是建符号数学工具箱是建立在功能强大的立在功能强大的 Maple 软件的基础上的，当软件的基础上的，当 Matlab 进进行符号运算时，它就请求行符号运算时，它就请求 Maple 软件去计算并将结果返软件去计算并将结果返回给回给 Mat

2、lab。qMatlab的符号数学工具箱可以完成几乎所有得符号运算功能。主要包括：符号表达式的运算，符号表达式的复合、化简，符号矩阵的运算，符号微积分、符号作图，符号代数方程求解，符号微分方程求解等。此外，该工具箱还支持可变精度运算，即支持以指定的精度返回结果。3Matlab符号运算特点 计算以计算以推理方式推理方式进行，因此不受计算误差累积所带来的进行，因此不受计算误差累积所带来的困扰。困扰。符号计算指令的调用比较简单，与数学教科书上的公式相近。符号计算可以给出完全正确的封闭解，或任意精度的数符号计算可以给出完全正确的封闭解，或任意精度的数值解（封闭解不存在时）。值解（封闭解不存在时）。符号计

3、算所需的运行时间相对较长。4Matlab符号运算举例 求一元二次方程求一元二次方程 ax2+bx+c=0 的根的根 solve(a*x2+b*x+c=0)求的根求的根 f(x)=(cosx)2 的一次导数的一次导数 x=sym(x);diff(cos(x)2)计算计算 f(x)=x2 在区间在区间 a,b 上的定积分上的定积分 syms a b x;int(x2,a,b)5q在进行符号运算时，必须先定义基本的符号对象，可以是符号常量、符号变量、符号表达式等。符号对象是一种数据结构。符号对象与符号表达式q含有符号对象的表达式称为符号表达式，Matlab在内部把符号表达式表示成字符串，以与数字变量

4、或运算相区别。q符号矩阵/数组：元素为符号表达式的矩阵/数组。6 sym 函数用来建立单个符号变量，一般调用格式为：q 符号对象的建立：符号对象的建立：sym 和和 syms符号对象的建立例：a=sym(a)符号变量=sym(A)参数A可以是一个数或数值矩阵，也可以是字符串a是符号变量b是符号常量 b=sym(1/3)C是符号矩阵 C=sym(1 ab;c d)7q 符号对象的建立：符号对象的建立：sym 和和 syms符号对象的建立 syms 命令用来建立多个符号变量，一般调用格式为：syms 符号变量符号变量1 符号变量符号变量2.符号变量符号变量n 例：syms a b c a=sym(

5、a);b=sym(b);c=sym(c);8q 符号表达式的建立：符号表达式的建立：例：建立符号表达式通常有以下2种方法：(1)用sym函数直接建立符号表达式。(2)使用已经定义的符号变量组成符号表达式。y=sym(sin(x)+cos(x)x=sym(x);y=sin(x)+cos(x)符号表达式的建立 syms x;y=sin(x)+cos(x)9Matlab符号运算采用的运算符和基本函数，在形状、名称和使用上，都与数值计算中的运算符和基本函数完全相同符号对象的基本运算q 基本运算符基本运算符普通运算：+-*/数组运算：.*./.矩阵转置：.例：X=sym(x11,x12;x21,x22;

6、x31,x32);Y=sym(y11,y12,y13;y21,y22,y23);Z1=X*Y;Z2=X.*Y;10符号对象的基本运算sin、cos、tan、cot、sec、csc、asin、acos、atan、acot、asec、acsc、exp、log、log2、log10、sqrtabs、conj、real、imagrank、det、inv、eig、lu、qr、svddiag、triu、tril、expm三角函数与反三角函数、指数函数、对数函数等q 基本函数基本函数11q 查找符号表达式中的符号变量查找符号表达式中的符号变量若表达式中有两个符号变量与x的距离相等，则ASCII码大者优先。查

7、找符号变量findsym(expr)按字母顺序列出符号表达式expr中的所有符号变量findsym(expr,N)按顺序列出expr中离x最近的N个符号变量常量pi,i,j不作为符号变量12例：f=sym(2*w-3*y+z2+5*a)findsym(f)findsym(f,3)findsym(f,1)findsym举例13符号表达式的替换subs(f,x,a)用用 a 替换字符函数替换字符函数 f 中的字符变量中的字符变量 x a是可以是是可以是 数数/数值变量数值变量/表达式表达式 或或 字符变量字符变量/表达式表达式若x是一个由多个字符变量组成的数组或矩阵，则a应该具有与x相同的形状的数

8、组或矩阵。q 用给定的用给定的数据数据替换符号表达式中的指定的替换符号表达式中的指定的符号变量符号变量14subs举例 f=sym(2*u);subs(f,u,2)f2=subs(f,u,u+2)a=3;subs(f2,u,a+2)subs(f2,u,a+2)syms x y f3=subs(f,u,x+y)subs(f3,x,y,1,2)ans=4f2=2*(u+2)ans=14ans=2*(a+2)+2)f3=2*x+2*yans=6 例：指出下面各条语句的输出结例：指出下面各条语句的输出结果果f=2*u15符号符号矩阵矩阵 A=sym(1+x,sin(x);5,exp(x)使用sym函数

9、直接生成 将数值矩阵转化成符号矩阵 符号矩阵中元素的引用和修改 B=2/3,sqrt(2);5.2,log(3);C=sym(B)A=sym(1+x,sin(x);5,exp(x);A(1,2)%引用引用 A(2,2)=sym(cos(x)%重新赋值重新赋值16六六类常见符号运算类常见符号运算q 因式分解、展开、合并、简化及通分等因式分解、展开、合并、简化及通分等q 计算极限计算极限q 计算导数计算导数q 计算积分计算积分q 符号求和符号求和q 代数方程和微分方程求解代数方程和微分方程求解17因式分解因式分解factor(f)syms x;f=x6+1;factor(f)lfactor 也可用

10、于正整数的分解 s=factor(100)factor(sym(12345678901234567890)l 大整数的分解要转化成符号常量大整数的分解要转化成符号常量18函数展开函数展开expand(f)syms x;f=(x+1)6;expand(f)l多项式展开l三角函数展开 syms x y;f=sin(x+y);expand(f)19合并同类项合并同类项collect(f,v):按指定变量按指定变量 v 进行合并进行合并collect(f):按按默认变量默认变量进行合并进行合并 syms x y;f=x2*y+y*x-x2+2*x;collect(f)collect(f,y)20函数简

11、化函数简化y=simple(f):对对 f 尝试多种不同的算法进行尝试多种不同的算法进行简化，返回其中最简短的形式简化，返回其中最简短的形式How,y=simple(f):y 为为 f 的最简短形式，的最简短形式，How 中记录的为简化过程中使用的方法。中记录的为简化过程中使用的方法。fRHOW2*cos(x)2-sin(x)23*cos(x)2-1 simplify(x+1)*x*(x-1)x3-xcombine(trig)x3+3*x2+3*x+1(x+1)3factorcos(3*acos(x)4*x3-3*xexpand21函数简化函数简化y=simplify(f):对对 f 进行简化

12、进行简化 syms x;f=sin(x)2+cos(x)2;simplify(f)syms c alpha beta;f=exp(c*log(sqrt(alpha+beta);simplify(f)22函数简化举例 syms x;f=(1/x3+6/x2+12/x+8)(1/3);y1=simplify(f)g1=simple(f)g2=simple(g1)l 多次使用多次使用 simple 可以达到最简表达。可以达到最简表达。例：简化23分式通分函数简化N,D=numden(f):N 为通分后的分子，为通分后的分子，D 为通分后的分母为通分后的分母 syms x y;f=x/y+y/x;N,

13、D=numden(f)n,d=numden(sym(112/1024)24horner多项式horner 多项式：嵌套形式的多项式 syms x;f=x4+2*x3+4*x2+x+1;g=horner(f)例：例：25计算极限limit(f,x,a):计算计算limit(f,a):当当默认变量默认变量趋向于趋向于 a 时的极限时的极限limit(f):计算计算 a=0 时的极限时的极限limit(f,x,a,right):计算右极限计算右极限limit(f,x,a,left):计算左极限计算左极限例：计算，syms x h n;L=limit(log(x+h)-log(x)/h,h,0)M=l

14、imit(1-x/n)n,n,inf)26计算导数g=diff(f,v)：求符号表达式求符号表达式 f 关于关于 v 的导数的导数g=diff(f)：求符号表达式求符号表达式 f 关于关于默认变量默认变量的导数的导数g=diff(f,v,n)：求求 f 关于关于 v 的的 n 阶导数阶导数q diff syms x;f=sin(x)+3*x2;g=diff(f,x)27计算积分int(f,v,a,b):计算定积分计算定积分int(f,a,b):计算关于计算关于默认变量默认变量的定积分的定积分int(f,v):计算不定积分计算不定积分int(f):计算关于计算关于默认变量默认变量的不定积分的不定

15、积分 syms x;f=(x2+1)/(x2-2*x+2)2;I=int(f,x)K=int(exp(-x2),x,0,inf)例：计算和28符号求和 syms n;f=1/n2;S=symsum(f,n,1,inf)S100=symsum(f,n,1,100)symsum(f,v,a,b):求和求和symsum(f,a,b):关于关于默认变量默认变量求和求和例：计算级数及其前100项的部分和例：计算函数级数 syms n x;f=x/n2;S=symsum(f,n,1,inf)29代数方程求解solve(f,v)：求方程关于指定自变量的解，求方程关于指定自变量的解，f 可以是可以是用字符串表

16、示的方程用字符串表示的方程、符号表达式符号表达式或或符号方程符号方程；l solve 也可解方程组也可解方程组(包含非线性包含非线性)；l 得不到解析解时，给出数值解。得不到解析解时，给出数值解。30微分方程微分方程求解求解q dsolvey=dsolve(eq1,eq2,.,cond1,cond2,.,v)其中其中 y 为输出的解，为输出的解，eq1、eq2、.为微分方程，为微分方程，cond1、cond2、.为初值条件，为初值条件，v 为自变量为自变量例例 1：求微分方程求微分方程 的通解，并验证。的通解，并验证。y=dsolve(Dy+2*x*y=x*exp(-x2),x)syms x;

17、diff(y)+2*x*y-x*exp(-x2)31微分方程微分方程求解求解q 几点说明几点说明l 如果省略初值条件，则表示求通解；如果省略初值条件，则表示求通解；l 如果省略自变量，则默认自变量为如果省略自变量，则默认自变量为 t dsolve(Dy=2*x,x);dy/dx=2xdsolve(Dy=2*x);dy/dt=2xl 若找不到解析解，则返回其积分形式。若找不到解析解，则返回其积分形式。l 微分方程中用微分方程中用 D 表示对表示对 自变量自变量 的导数，如：的导数，如：Dy y；D2y y；D3y y32微分方程微分方程求解求解例例 2：求微分方程求微分方程 满足初值条件满足初值

18、条件 的特解，并画出解函数的图形。的特解，并画出解函数的图形。y=dsolve(x*Dy+y-exp(x)=0,.y(1)=2*exp(1),x)ezplot(y);33微分方程微分方程求解求解例例3：求微分方程组求微分方程组 在初值条件在初值条件 下的特解，并画出解函数的图形。下的特解，并画出解函数的图形。x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0,.x(0)=1,y(0)=0,t)ezplot(x,y,0,1.3);注：解微分方程组时，如果所给的输出个数与方程个数相同，注：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解则方程组的解 按词典顺序

19、按词典顺序 输出。输出。34其它运算其它运算反函数finverse(f,v)：求求 f 关于指定变量关于指定变量 v 的反函数的反函数finverse(f)：求求 f 关于默认变量的反函数关于默认变量的反函数 syms x t;f=x2+2*t;g1=finverse(f,x)g2=finverse(f,t)例：计算函数的反函数35MATLAB应用应用-解线性规划解线性规划36用用MATLAB优化工具箱解线性规划优化工具箱解线性规划minz=cX 1、模型：命令：x=linprog（c，A，b）2、模型：minz=cX 命令：x=linprog（c，A，b，Aeq,beq）注意：若没有不等式：

20、存在，则令A=，b=.373、模型：minz=cX VLBXVUB命令：1x=linprog（c，A，b，Aeq,beq,VLB，VUB）2x=linprog（c，A，b，Aeq,beq,VLB，VUB,X0）注意：1若没有等式约束:,则令Aeq=,beq=.2其中X0表示初始点4、命令：x,fval=linprog()返回最优解及处的目标函数值fval.38解解 编写编写M文件文件xxgh1.m如下：如下：c=-0.4-0.28-0.32-0.72-0.64-0.6;A=0.01 0.01 0.01 0.03 0.03 0.03;0.02 0 0 0.05 0 0;0 0.02 0 0 0.

21、05 0;0 0 0.03 0 0 0.08;b=850;700;100;900;Aeq=;beq=;vlb=0;0;0;0;0;0;vub=;x,fval=linprog(c,A,b,Aeq,beq,vlb,vub)39解解:编写编写M文件文件xxgh2.m如下：如下：c=6 3 4;A=0 1 0;b=50;Aeq=1 1 1;beq=120;vlb=30,0,20;vub=;x,fval=linprog(c,A,b,Aeq,beq,vlb,vub)40MATLAB应用应用 解非线性规划解非线性规划41一、二次规划一、二次规划(Quadratic Program)概念概念42二、二、Mat

22、lab中求解二次规划中求解二次规划4344转转化化为为matlabmatlab求解格式：求解格式：4546 定义定义如果目标函数或约束条件中至少有一个是非线性函数,则最优化问题就叫做非线性规划问题非线性规划问题四、非线性规划的基本概念四、非线性规划的基本概念 一般形式一般形式:（1）其中 ，是定义在 R Rn 上的实值函数()nTnRxxxX=,21L()()=.,.,2,10m;1,2,.,0.ljXhiXgtsji47五、非线性规划的基本解法五、非线性规划的基本解法SUTM外点法外点法SUTM内点法（障碍罚函数法）内点法（障碍罚函数法）1 罚函数法罚函数法2 近似线性规划法近似线性规划法4

23、8 1 1、罚函数法罚函数法 罚函数法罚函数法基本思想是通过构造罚函数把约束问题转化为一系列无约束最优化问题，进而用无约束最优化方法去求解这类方法称为序列无约束最小化方法序列无约束最小化方法简称为SUMTSUMT法法 其一为SUMTSUMT外点法外点法，其二为SUMTSUMT内点法内点法49Matlab求解非线性规划问题求解非线性规划问题 其中其中X为为n维变元向量，维变元向量，G(X)与与Ceq(X)均为非线性函数组成的均为非线性函数组成的向量。向量。50 1 首先建立首先建立M文件文件fun.m,用来定义目标函数用来定义目标函数F（X）:function f=fun(X);f=F(X);M

24、ATLAB求解上述问题，基本步骤分三步求解上述问题，基本步骤分三步513 建立主程序建立主程序.求解非线性规划的函数是求解非线性规划的函数是fmincon,命令的基本命令的基本格式如下：格式如下：(1)x=fmincon(fun,X0,A,b)(2)x=fmincon(fun,X0,A,b,Aeq,beq)(3)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB)(4)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon)(5)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon,options)(6

25、)x,fval=fmincon()(7)x,fval,exitflag=fmincon()(8)x,fval,exitflag,output=fmincon()输出极值点M文件迭代的初值参数说明变量上下限fmincon函数可能会给出局部最优解，这与初值函数可能会给出局部最优解，这与初值X0的选取有关的选取有关521写成标准形式写成标准形式：s.t.2x1+3x2 6 s.t.x1+4x2 5 x1,x2 0例例532先建立先建立M-文件文件 fun3m:function f=fun3(x);f=-x(1)-2*x(2)+(1/2)*x(1)2+(1/2)*x(2)23再建立主程序youh2m：

26、x0=1;1;A=2 3;1 4;b=6;5;Aeq=;beq=;VLB=0;0;VUB=;x,fval=fmincon(fun3,x0,A,b,Aeq,beq,VLB,VUB)4运算结果为：运算结果为：x=07647 10588 fval=-20294541 1先建立先建立M文件文件fun4m定义目标函数定义目标函数:function f=fun4(x);f=exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1);x1+x2=0s.t.1.5+x1x2-x1-x20-x1x210 0例例 2再建立再建立M文件文件myconm定义非线性约束：定义非线性约

27、束：function g,ceq=mycon(x)g=x(1)+x(2);15+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10;553主程序主程序youh3m为为:x0=-1;1;A=;b=;Aeq=1 1;beq=0;vlb=;vub=;x,fval=fmincon(fun4,x0,A,b,Aeq,beq,vlb,vub,mycon)4运算结果为运算结果为：x=-12250 12250 fval=1895156 例例1先建立先建立M文件文件funm定义目标函数定义目标函数:function f=fun(x);f=-2*x(1)-x(2);2再建立再建立M文件文件mycon

28、2m定义非线性约束：定义非线性约束：function g,ceq=mycon2(x)g=x(1)2+x(2)2-25;x(1)2-x(2)2-7;573 主程序主程序fxxm为为:x0=3;25;VLB=0 0;VUB=5 10;x,fval,exitflag,output =fmincon(fun,x0,VLB,VUB,mycon2)58MATLAB应用应用 微分方程解法微分方程解法59要求目的要求目的主要内容主要内容2、学会用、学会用Matlab求微分方程的数值解求微分方程的数值解.1、学会用、学会用Matlab求简单微分方程的解析解求简单微分方程的解析解.1 1、求简单微分方程的解析解求

29、简单微分方程的解析解.2、求微分方程的数值解、求微分方程的数值解.60求微分方程的数值解求微分方程的数值解（一）常微分方程数值解的定义一）常微分方程数值解的定义（二）建立数值解法的一些途径（二）建立数值解法的一些途径（三）用三）用Matlab软件求常微分方程的数值解软件求常微分方程的数值解返回61微分方程的解析解微分方程的解析解求微分方程（组）的解析解命令:dsolve(方程方程1,方程方程2,方程方程n,初始条件初始条件,自变量自变量)结果：u=tg(t-c)62 解解输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)结果为:y=3e-2xsin

30、（5x）63解解输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z,t)；x=simple(x)%将x化简y=simple(y)z=simple(z)结果为：x=(c1-c2+c3+c2e-3t-c3e-3t)e2ty=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t64微分方程的数值解微分方程的数值解（一）常微分方程数值解的定义一）常微分方程数值解的定义在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解。而在实际上对初

31、值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。因此，研究常微分方程的数值解法是十分必要的因此，研究常微分方程的数值解法是十分必要的。65（二）建立数值解法的一些途径二）建立数值解法的一些途径1、用差商代替导数、用差商代替导数若步长h较小，则有故有公式：此即欧拉法欧拉法。662、使用数值积分、使用数值积分对方程y=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：实际应用时，与欧拉公式结合使用：此即改进的欧拉法改进的欧拉法。故有公式：673、使用泰勒公式、使用泰勒公式以此方法为基础，有龙格龙格-库塔法库塔法、线性多步法线性多步

32、法等方法。4、数值公式的精度、数值公式的精度当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。68（三）用三）用Matlab软件求常微分方程的数值解软件求常微分方程的数值解t，x=solver（f,ts,x0,options）ode45ode23ode113ode15sode23s由待解方程写成的m-文件名ts=t0，tf，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库

33、塔-芬尔格算法ode45：运用组合的4/5阶龙格-库塔-芬尔格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（reltol,rt,abstol,at）,rt，at：分别为设定的相对误差和绝对误差.691、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意注意:70解解:令 y1=x，y2=y11、建立m-文件vdp1000.m如下：function dy=vdp1000(t,y)dy=zeros(

34、2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)2)*y(2)-y(1);2、取t0=0，tf=3000，输入命令：T,Y=ode15s(vdp1000,03000,20);plot(T,Y(:,1),-)3、结果如图71解解 1、建立m-文件rigid.m如下：function dy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);2、取t0=0，tf=12，输入命令：T,Y=ode45(rigid,0 12,0 1 1);plot(T,Y(:,1),-,T,Y(

35、:,2),*,T,Y(:,3),+)3、结果如图图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.72Matlab的应用的应用 -插值和拟合插值和拟合73拉格朗日插值拉格朗日插值分段线性插值分段线性插值三次样条插值三次样条插值一一 维维 插插 值值一、一、插值的定义插值的定义二、插值的方法二、插值的方法三、用三、用Matlab解插值问题解插值问题74二维插值二维插值一、一、二维插值定义二维插值定义二、二、网格节点插值法网格节点插值法三、三、用用MatlabMatlab解插值问题解插值问题最邻近插值最邻近插值分片线性插值分片线性插值双线性插值双线性插值网格节点数据的插值网格节点

36、数据的插值散点数据的插值散点数据的插值75一维插值的定义一维插值的定义已知已知 n+1个节点个节点其中其中互不相同，不妨设互不相同，不妨设求任一插值点求任一插值点处的插值处的插值节点可视为由节点可视为由产生产生,，表达式复杂表达式复杂,，或无封闭形式或无封闭形式或未知或未知76 构造一个构造一个(相对简单的相对简单的)函数函数通过全部节点通过全部节点,即即再用再用计算插值，即计算插值，即77称为拉格朗日插值基函数拉格朗日插值基函数。已知函数f(x)在n+1个点x0,x1,xn处的函数值为y0,y1,yn。求一n次多项式函数Pn(x)，使其满足：Pn(xi)=yi,i=0,1,n.解决此问题的拉

37、格朗日插值多项式公式如下其中Li(x)为n次多项式：拉格朗日拉格朗日(Lagrange)插值插值78拉格朗日拉格朗日(Lagrange)插值插值特别地特别地:两点一次两点一次(线性线性)插值多项式插值多项式:三点二次三点二次(抛物抛物)插值多项式插值多项式:79分段线性插值分段线性插值计算量与n无关;n越大，误差越小.xjxj-1xj+1x0 xnxoy80比分段线性插值更光滑。比分段线性插值更光滑。xyxi-1 xiab在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次

38、样条插值就是一个很好的例子。三次样条插值三次样条插值81 三次样条插值82用用MATLABMATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插最邻近插值值linear ：线性插线性插值；值；spline ：三次样条三次样条插值；插值；cubic ：立方插值。立方插值。缺省时：缺省时：分段线性插值。分段线性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范

39、围。83 例：在例：在1-121-12的的1111小时内，每隔小时内，每隔1 1小时测量一次温小时测量一次温度，测得的温度依次为：度，测得的温度依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424。试估计每隔。试估计每隔1/101/10小时的温小时的温度值。度值。hours=1:12;temps=589152529313022252724;h=1:0.1:12;t=interp1(hours,temps,h,spline);(直接输出数据将是很多的)plot(hours,temps,+,h,t,hours,temps,r:)%

40、作图xlabel(Hour),ylabel(DegreesCelsius)84二维插值的定义二维插值的定义 xyO O第一种（网格节点）：第一种（网格节点）：85 已知已知 m n个节点个节点 其中其中互不相同，不妨设互不相同，不妨设 构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，即86第二种（散乱节点）：第二种（散乱节点）：yx0 087已知已知n个节点个节点其中其中互不相同，互不相同，构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，即返回返回88 注意：注意：最邻近插值一般不连续。具有

41、连续性的最简单的插值是分片线性插值。最邻近插值最邻近插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。89将四个插值点（矩形的四个顶点）处的函数值依次简记为：分片线性插值分片线性插值xy (xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)O Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f490插值函数为：第二片(上三角形区域)：(x,y)满足插值函数为：注意注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分

42、片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足91双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值双线性插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O92 要求要求x0,y0 x0,y0单调；单调；x x，y y可取可取为矩阵，或为矩阵，或x x取取行向量，行向量，y y取为列向量，取为列向量，x,yx,y的值分别不能超出的值分别不能超出x0,y0 x0,y0的范围。的范围。z=interp2(x0

43、,y0,z0,x,y,method)被插值点插值方法用用MATLAB作网格节点数据的插值作网格节点数据的插值插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值缺省时缺省时,双线性插值双线性插值93例：测得平板表面例：测得平板表面3*53*5网格点处的温度分别为：网格点处的温度分别为：82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 84 8484 82 85 86 82 85 86 试作出平板表面的温度

44、分布曲面试作出平板表面的温度分布曲面z=f(x,y)z=f(x,y)的图形。的图形。输入以下命令：x=1:5;y=1:3;temps=8281808284;7963616581;8484828586;mesh(x,y,temps)1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.2以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.94再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图.95通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果

45、进行比较。96 插值函数插值函数griddata格式为格式为:cz=griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作散点数据的插值计算作散点数据的插值计算 要求要求cxcx取行向量，取行向量，cycy取为列向量取为列向量。被插值点插值方法插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值v4-Matlab提供的插值方法提供的插值方法缺省时缺省时,双线性插值双线性插值97拟拟 合合2.2.拟合的基本原理拟合的基本原理1.拟合问题引例拟合问题引

46、例98拟拟 合合 问问 题题 引引 例例 1 1温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032已知热敏电阻数据：已知热敏电阻数据：求求60600C时的电阻时的电阻R。设设R=at+ba,b为待定系数为待定系数99拟拟 合合 问问 题题 引引 例例 2 2t(h)0.25 0.5 1 1.5 2 3 4 6 8c(g/ml)19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注射300m

47、g)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).作半对数坐标系作半对数坐标系(semilogy)下的图形下的图形100曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi)i=1,n,寻求一个函数（曲线）寻求一个函数（曲线）y=f(x),使使 f(x)在在某种准则下与所有某种准则下与所有数据点最为接近，即曲线拟合得最好。数据点最为接近，即曲线拟合得最好。+xyy=f(x)(xi,yi)i i 为点为点（xi,yi)与与曲线曲线 y=f(x)的的距离距离101拟合与插值的关系拟合与插值的关系函数

48、插值与曲线拟合都是要根据一组数据构造一个函数作为近函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上完全不同。似，由于近似的要求不同，二者的数学方法上完全不同。实例：实例：下面数据是某次实验所得，希望得到x和f之间的关系？问题：问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，就是数据拟合数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所有数据点，就是插值问题插值问题；102最临近插值、线性插值、样条插值与曲线拟合结果：最临近插值、线性插值、样条插

49、值与曲线拟合结果：103曲线拟合问题最常用的解法曲线拟合问题最常用的解法线性最小二乘法的基本思路线性最小二乘法的基本思路第一步:先选定一组函数先选定一组函数r1(x),r2(x),rm(x),mn,令令f(x)=a1r1(x)+a2r2(x)+amrm(x)（1）其中其中a1,a2,am为待定系数。为待定系数。第二步:确定确定a1,a2,am的的准则（最小二乘准则）：准则（最小二乘准则）：使使n个点个点（xi,yi)与与曲线曲线 y=f(x)的的距离距离 i 的平方和最小的平方和最小。记记问题归结为，求问题归结为，求a1,a2,am 使使 J(a1,a2,am)最小。最小。104线性最小二乘法

50、的求解：预备知识线性最小二乘法的求解：预备知识超定方程组超定方程组：方程个数大于未知量个数的方程组：方程个数大于未知量个数的方程组即即 Ra=y其中其中超定方程一般是不存在解的矛盾方程组。超定方程一般是不存在解的矛盾方程组。如果有向量如果有向量a使得使得 达到最小，达到最小，则称则称a为上述为上述超定方程的最小二乘解超定方程的最小二乘解。105线性最小二乘法的求解线性最小二乘法的求解定理：定理：当当R RT TR R可逆时，超定方程组（可逆时，超定方程组（3 3）存在最小二乘解，）存在最小二乘解，且即为方程组且即为方程组 R RT TRaRa=R RT Ty y -正则（正规）方程组正则（正规