第3章差分微分matlab.ppt

微分方程与差分方程上机合肥师范学院 MATLAB的符号运算的符号运算 符号计算是对未赋值的符号对象(可以是常数、变量、表达式)进行运算和处理。MATLAB具有符号数学工具箱(Symbolic Math Toolbox)，将符号运算结合到MATLAB的数值运算环境。符号数学工具箱是建立在Maple软件基础上的。符号表达式的建立符号表达式的建立创创建符号常量建符号常量符号常量是不含变量的符号表达式，用sym命令来创建符号常量。语法：语法：sym(常量常量)%创建符号常量例如：a=sym(sin(2)a=sin(2)创建符号变量和表达式创建符号变量和表达式1.使用使用sym命令创建符号变量和表达式命令创建符号变量和表达式语法：语法：sym(表达式表达式)%创建符号表达式 符号变量名符号变量名=sym(表达式表达式)%符号表达式赋给 符号变量2.使用使用syms命令创建符号变量和符号表达式命令创建符号变量和符号表达式syms用于创建多个多个符号变量语法：语法：syms(arg1,arg2,参数参数)%把字符变量定义为符号变量syms arg1 arg2,参数参数%把字符变量定义为符号变量的简洁形式简洁形式【例例】使用syms命令创建符号变量和符号表达式。syms a b c x%创建多个符号变量创建多个符号变量 f2=a*x2+b*x+c%创建符号表达式创建符号表达式 f2=a*x2+b*x+c syms(a,b,c,x)f3=a*x2+b*x+c;%创建符号表达式创建符号表达式 符号矩阵符号矩阵用sym和syms命令也可以创建符号矩阵。例如，A=sym(a,b;c,d)A=a,b c,d syms a b c d A=a b;c d A=a,b c,d 符号表达式的代数运算符号表达式的代数运算符号运算与数值运算的区别主要有以下几点：符号运算与数值运算的区别主要有以下几点：符号运算与数值运算的区别主要有以下几点：符号运算与数值运算的区别主要有以下几点：符号运算不需要进行数值运算，不会出现截断误差，因此符号运算是非常准确的。符号运算可以得出完全的封闭解或任意精度的数值解。符号运算的时间较长，而数值型运算速度快。符号表达式的代数运算符号表达式的代数运算1.符号运算中的运算符符号运算中的运算符(1)基本运算符运算符“”，“”，“*”，“”，“/”，“”分别实现符号矩阵的加、减、乘、左除、右除、求幂运算。运算符“.*”，“./”，“.”，“.”分别实现符号数组的乘、除、求幂，即数组间元素与元素的运算。运算符“”，“.”分别实现符号矩阵的共轭转置、非共轭转置。(2)关系运算符在符号对象的比较中，没有没有“大于”、“大于等于”、“小于”、“小于等于”的概念，而只有是否“等于”的概念。运算符“=”、“=”分别对运算符两边的符号对象进行“相等”、“不等”的比较。当为“真”时，比较结果用1表示；当为“假”时，比较结果则用0表示。2.函数运算函数运算(1)三角函数和双曲函数三角函数包括sin、cos、tan；双曲函数包括sinh、cosh、tanh；三角反函数除了atan2函数仅能用于数值计算外，其余的asin、acos、atan函数在符号运算中与数值计算的使用方法相同。(2)指数和对数函数指数函数sqrt、exp的使用方法与数值计算的完全相同；对数函数在符号计算中只有自然对数log(表示ln)，而没有数值计算中的log2和log10。(3)复数函数复数的共轭conj、求实部real、求虚部imag和求模abs函数与数值计算中的使用方法相同。但注意，在符号计算中，MATLAB没有提供求相角的命令。(4)矩阵代数命令MATLAB提供的常用矩阵代数命令有diag，triu，tril，inv，det，rank，poly，eig、expm等，它们的用法几乎与数值计算中的情况完全一样。【例例】求矩阵的行列式值、非共轭转置和特征值。syms a11 a12 a21 a22 A=a11 a12;a21 a22%创建符号矩阵创建符号矩阵 A=a11,a12 a21,a22 det(A)%计算行列式计算行列式 ans=a11*a22-a12*a21 A.%计算非共轭转置计算非共轭转置 ans=a11,a21 a12,a22 eig(A)%计算特征值计算特征值 【例例】符号表达式f=2x2+3x+4与g=5x+6的代数运算。f=sym(2*x2+3*x+4)f=2*x2+3*x+4 g=sym(5*x+6)g=5*x+6 f+g%符号表达式相加符号表达式相加 ans=2*x2+8*x+10 f*g%符号表达式相乘符号表达式相乘 ans=(2*x2+3*x+4)*(5*x+6)1.自由变量的确定原则自由变量的确定原则小写字母i和j不能作为自由变量。符号表达式中如果有多个字符变量，则按照以下顺序选择自由变量：首先选择x作为自由变量；如果没有x，则选择在字母顺序中最接近x的字符变量；如果与x相同距离，则在x后面的优先。大写字母比所有的小写字母都靠后。2.findsym函数函数如果不确定符号表达式中的自由符号变量，可以用findsym函数来自动确定。语法：语法：findsym(f,n)%确定自由符号变量 说明：f可以是符号表达式或符号矩阵；n为按顺序得出符号变量的个数，当n省略时，则不按顺序得出f中所有的符号变量。符号表达式的操作和转换符号表达式的操作和转换1、符号表达式中自由变量的确定、符号表达式中自由变量的确定 2、符号表达式的化简、符号表达式的化简(1)pretty函数函数 将给出排版形式的输出结果。(2)collect函数函数 将表达式中相同次幂的项合并，也可以将表达式中相同次幂的项合并，也可以再输入一个参数指定以哪个变量的幂次合并。再输入一个参数指定以哪个变量的幂次合并。(3)expand函数函数 将表达式展开成多项式形式。将表达式展开成多项式形式。(4)horner函数函数 将表达式转换为嵌套格式。将表达式转换为嵌套格式。(5)factor函数函数 将表达式转换为嵌套格式。将表达式转换为嵌套格式。(6)simplify函数函数 利用函数规则对表达式进行化简。利用函数规则对表达式进行化简。(7)simple函数函数 调用MATLAB的其他函数对表达式进行综合化简，并显示化简过程。3、符号表达式的替换、符号表达式的替换MATLAB 中，可以通过符号替换使表达式的形式简化。符号工具箱中提供了两个函数用于表达式的替换：1subexpr该函数自动将表达式中重复出现的比较长的子表达式或字符串用变量替换，该函数的调用格式为：subexpr(s,s1)，指定用符号变量 s1 来代替符号表达式s（可以是矩阵）中重复出现的字符串。替换后的结果由 ans 返回，被替换的字符串由 s1返回；Y,s1=subexpr(X,s1)，该命令与上面的命令不同之处在于第二个参数为字符串，该命令用来替换表达式中重复出现的字符串。2.subs函数 subs 可以用指定符号替换表达式中的某一特定符号。subs(s)subs(s,new)subs(s,old,new)4、求反函数和复合函数、求反函数和复合函数语法：语法：finverse(f,v)%对指定自变量v的函数f(v)求反函数 compose(f,g)%计算复合函数f(g(x)5、符号表达式与多项式的转符号表达式与多项式的转换换 构成多项式的符号表达式f(x)可以与多项式系数构成的行向量进行相互转换，MATLAB提供了函数sym2poly和poly2sym实现相互转换。【例例1】将符号表达式2x+3x2+1转换为行向量。f=sym(2*x+3*x2+1)f=2*x+3*x2+1 sym2poly(f)%转换为按降幂排列的行向量转换为按降幂排列的行向量ans=3 2 1 【例例2】将行向量转换为符号表达式。g=poly2sym(1 3 2)%默认默认x为符号变量的符号表达式为符号变量的符号表达式g=x2+3*x+2 符号极限、微积分和级数求和符号极限、微积分和级数求和符符号极限号极限 假定符号表达式的极限存在，Symbolic Math Toolbox提供了直接求表达式极限的函数limit，函数limit的基本用法如表所示。符号微分符号微分函数diff是用来求符号表达式的微分。语法：语法：diff(f)%求f对自由变量的一阶微分 diff(f,t)%求f对符号变量t的一阶微分 diff(f,n)%求f对自由变量的n阶微分 diff(f,t,n)%求f对符号变量t的n阶微分符号积分符号积分 积分有定积分和不定积分，运用函数int可以求得符号表达式的积分。语法：语法：int(f,t)%求符号变量t的不定积分int(f,t,a,b)%求符号变量t的积分int(f,t,m,n)%求符号变量t的积分说明说明：t为符号变量，当t省略则为默认自由变量；a和b为数值，a,b为积分区间；m和n为符号对象，m,n为积分区间；与符号微分相比，符号积分复杂得多。因为函数的积分有时可能不存在，即使存在，也可能限于很多条件，MATLAB无法顺利得出。当MATLAB不能找到积分时，它将给出警告提示并返回该函数的原表达式。符号级数符号级数1.symsum函数函数语法：语法：symsum(s,x,a,b)%计算表达式s的级数和.说明：x为自变量，x省略则默认为对自由变量求和；s为符号表达式；a,b为参数x的取值范围。2.taylor函数函数语法：语法：taylor(F,x,n)%求泰勒级数展开说明：x为自变量，F为符号表达式；对F进行泰勒级数展开至n项，参数n省略则默认展开前5项。【例例】求级数 1+x+x2+xk+的和。syms x k s1=symsum(1/k2,1,10)%计算级数的前计算级数的前10项和项和 s1=1968329/1270080 s2=symsum(1/k2,1,inf)%计算级数和计算级数和 s2=1/6*pi2 s3=symsum(xk,k,0,inf)%计算对计算对k为自变为自变量的级数和量的级数和 s3=-1/(x-1)【例例】求ex的泰勒展开式 syms xs1=taylor(exp(x),8)%展开前展开前8项项 s1=1+x+1/2*x2+1/6*x3+1/24*x4+1/120*x5+1/720*x6+1/5040*x7 s2=taylor(exp(x)%默认展开前默认展开前5项项 s2=1+x+1/2*x2+1/6*x3+1/24*x4+1/120*x5 符号方程的求解符号方程的求解代代数方程数方程语法：语法：solve(eq,v)%求方程关于指定变量的解solve(eq1,eq2,v1,v2,)%求方程组关于指定变量的解 说说明明：eq可以是含等号的符号表达式的方程，也可以是不含等号的符号表达式，但所指的仍是令eq=0的方程；当参数v省略时，默认为方程中的自由变量；其输出结果为结构数组类型。符号常微分方程符号常微分方程语法：语法：dsolve(eq,con,v)%求解微分方程dsolve(eq1,eq2,con1,con2,v1,v2)%求解微分方程组说明说明：eq为微分方程；con是微分初始条件，可省；v为指定自由变量，省略时则默认为x或t为自由变量；输出结果为结构数组类型。当y是因变量时，微分方程eq的表述规定为：y的一阶导数 表示为Dyy的n阶导数 表示为Dny微分初始条件con应写成y(a)=b，Dy(c)=d的格式。一、利用一、利用MatlabMatlab求微分方程的解析解求微分方程的解析解 求微分方程（组）的解析解命令:dsolve(方程方程1,方程方程2,方程方程n,初始条件初始条件,自变量自变量)结 果：u=tg(t-c)解解 输入命令:y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)结 果 为:y=3e-2xsin（5x）解解 输入命令：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)e2t y=-c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z=(-c1e-4t+c2e-4t+c1-c2+c3)e2t 建立数值解法的一些途径建立数值解法的一些途径1、用差商代替导数、用差商代替导数 若步长h较小，则有故有公式：此即欧拉法欧拉法。2、使用数值积分、使用数值积分对方程y=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：实际应用时，与欧拉公式结合使用：此即改进的欧拉法改进的欧拉法。故有公式：3、使用泰勒公式、使用泰勒公式 以此方法为基础，有龙格龙格-库塔（库塔（Runge Kutta）法）法、线性多步法线性多步法等方法。4、数值公式的精度、数值公式的精度 当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。返 回（三）可以用三）可以用Matlab软件求常微分方程的数值解软件求常微分方程的数值解t，x=solver（f,ts,x0,options）ode45 ode23 ode113ode15sode23s由待解方程写成的m-文件名ts=t0，tf，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔-芬尔格算法ode45：运用组合的4/5阶龙格-库塔-芬尔格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（reltol,rt,abstol,at）,rt，at：分别为设定的相对误差和绝对误差.不同求解器Solver的特点求解器SolverODE类型特点说明ode45非刚性一步算法；4，5阶Runge-Kutta方程；累计截断误差达(x)3大部分场合的首选算法ode23非刚性一步算法；2，3阶Runge-Kutta方程；累计截断误差达(x)3使用于精度较低的情形求解器SolverODE类型特点说明ode113非刚性多步法；Adams算法；高低精度均可到10-310-6计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gears反向数值微分；精度中等若ode45失效时，可尝试使用不同求解器Solver的特点求解器SolverODE类型特点说明ode23s刚性一步法；2阶Rosebrock算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短参数设置属性名取值含义AbsTol有效值：正实数或向量缺省值：1e-6绝对误差对应于解向量中的所有元素；向量则分别对应于解向量中的每一分量RelTol有效值：正实数缺省值：1e-3相对误差对应于解向量中的所有元素。在每步(第k步)计算过程中，误差估计为：e(k)=max(RelTol*abs(y(k),AbsTol(k)参数设置属性名取值含义NormControl有效值：on、off缺省值：off为on时，控制解向量范数的相对误差，使每步计算中，满足：norm(e)1缺省值：k=1若k1，则增加每个积分步中的数据点记录，使解曲线更加的光滑参数设置属性名取值含义Jacobian有效值：on、off缺省值：off若为on时，返回相应的ode函数的Jacobi矩阵Jpattern有效值：on、off缺省值：off为on时，返回相应的ode函数的稀疏Jacobi矩阵参数设置属性名取值含义Mass有效值：none、M、M(t)、M(t,y)缺省值：noneM：不随时间变化的常数矩阵M(t)：随时间变化的矩阵M(t,y)：随时间、地点变化的矩阵MaxStep有效值：正实数缺省值：tspans/10最大积分步长创建函数创建函数function2,保存在保存在function2.m中中function f=function2(t,x)f=-x.2;在命令窗口中执行在命令窗口中执行 t,x=ode45(function2,0,1,1);plot(t,x,-,t,x,o);xlabel(time t0=0,tt=1);ylabel(x values x(0)=1);创建函数创建函数function3,保存在保存在function3.m中：中：function f=function3(t,x)f=x(1)-0.1*x(1)*x(2)+0.01*t;.-x(2)+0.02*x(1)*x(2)+0.04*t;运行命令文件runf3.mt,x=ode45(function3,0,20,30;20);plot(t,x);xlabel(time t0=0,tt=20);ylabel(x values x1(0)=30,x2(0)=20);创建函数function4，存在function4.m中function f=function4(t,x)global Uf=x(2);U*(1-x(1)2)*x(2)-x(1);再运行命令文件再运行命令文件runf4.m：global UU=7;Y0=1;0;t,x=ode45(function4,0,40,Y0);x1=x(:,1);x2=x(:,2);plot(t,x1,t,x2)