数学软件 Matlab —— Matlab 符号运算
主要内容 Matlab 符号运算介绍 符号对象与基本符号运算 symvar、 subs 和 vpa 常见的符号计算(重点内容)
符号运算 符号运算的特点 计算以推理方式进行,不受计算误差累积所带来的困扰 符号计算可以给出完全正确的封闭解,或任意精度的数值解(封闭解不存在时 ) 符号计算指令的调用比较简单,与教科书上的公式相近 符号计算所需的运行时间相对较长
Matlab 符号运算 Matlab 符号运算 Matlab 符号运算是通过符号数学工具箱(Symbolic Math Toolbox)来实现的。 Matlab 的符号数学工具箱可以完成几乎所有得符号运算功能,如:符号表达式的运算,符号矩阵的运算,符号微积分,符号作图,符号代数方程求解,符号微分方程求解等。 此外,该工具箱还支持可变精度运算,即支持以指定的精度返回结果。
符号运算举例 solve('a*x^2+b*x+c=0') x=sym('x'); diff(cos(x)^2) syms a b x; 求的根 f (x) = (cos x)2 的一次导数 x=sym('x'); diff(cos(x)^2) 计算 f (x) = x2 在区间 [a, b] 上的定积分 syms a b x; int(x^2,a,b)
内容提要 Matlab 符号运算介绍 符号对象与基本符号运算 symvar、 subs 和 vpa 常见的符号计算
符号对象 Matlab 符号对象 在进行符号运算时,必须先定义基本的符号对象,可以是 符号变量、符号表达式等 符号对象是一种数据结构 在进行符号运算时,必须先定义基本的符号对象,可以是 符号变量、符号表达式等 符号对象是一种数据结构 符号表达式:含有符号对象的表达式称 符号矩阵/数组:元素为符号表达式的矩阵/数组
符号对象的建立 符号对象的定义/声明:sym、syms a 是符号变量 b 是符号常量 C 是符号矩阵 例: 符号变量 = sym(x) 参数 x 可以是一个数或数值矩阵,也可以是字符串 例: a=sym('a') a 是符号变量 b=sym('1/3') b 是符号常量 C=sym('[1 ab; c d]') C 是符号矩阵
符号对象的建立 符号对象的定义/声明:sym、syms syms 符号变量1 符号变量2 ... 符号变量n 例: a=sym('a'); b=sym('b'); c=sym('c'); 例: syms a b c;
符号表达式 符号表达式:含符号对象的表达式 例: 建立符号表达式通常有以下 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)
Matlab 符号运算采用的运算符和基本函数,在形状、名称和使用上,都与数值计算中的运算符和基本函数完全相同 基本符号运算 Matlab 符号运算采用的运算符和基本函数,在形状、名称和使用上,都与数值计算中的运算符和基本函数完全相同 基本运算 普通运算: + - * \ / ^ 数组运算: .* .\ ./ .^ 矩阵转置: ' .' 基本数学函数 三角函数与反三角函数、指数函数、对数函数等 sin,cos,asin,acos,exp,log,abs,diag,tril,triu, ...
符号矩阵 符号矩阵的生成 符号矩阵中元素的引用和修改 使用 sym 函数直接生成 A=sym('[1+x, sin(x); 5, exp(x)]') 将数值矩阵转化成符号矩阵 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)') % 重新赋值
内容提要 Matlab 符号运算介绍 符号对象与基本符号运算 symvar、 subs 和 vpa 常见的符号计算
symvar 列出符号表达式中的符号变量 symvar(s) symvar(s, N) 例: 按字母顺序列出符号表达式 s 中的所有符号变量 列出符号表达式 s 中离 x 最近的 N 个符号变量 若有两个符号变量与 x 的距离相等,则ASCII 码大者优先 常量 pi, i, j 不作为符号变量 例: f=sym('2*v-3*y+z^2+5*a') symvar(f) symvar(f,2)
subs 符号替换 subs(s,x,a) 例: 用给定的数据替换符号表达式中的指定的符号变量 用 a 替换符号表达式 s 中的符号变量 x Matlab 演示 syms x y u v; f1=2*x+y-1; f2=subs(f1,x,u) f3=subs(f1,y,2+3) f3=subs(f1,{x,y},{u,v})
subs 举例 例:指出下面各条语句的输出结果 f=sym('2*u'); f1=subs(f,'u',2) f2=subs(f,'u','u+2') f3=subs(f,'u',[1,2]) a=3; f4=subs(f2,'u',a+2) f5=subs(f2,'u','a+2') syms x y; f6=subs(f,'u',x+y) f7=subs(f6,{x,y},{1,2}) f8=subs(f6,{x,y},{x+y,x+y}) f=2*u f1=4 f2=2*u+4 f3=[2,4] f4=14 f5=2*a+8 f6=2*x+2*y f7=6 f8=4*x+4*y Matlab 演示 符号变量若没有声明,则需要加上单引号!
vpa 可变精度 vpa(s,n) 例: 计算表达式 s 的值,保留 n 位有效数字 返回值是符号对象 x1=vpa(sin(pi/2),10) x2=vpa(pi^3,3) x3=vpa(pi,100) Matlab 演示
内容提要 Matlab 符号运算介绍 符号对象与基本符号运算 symvar、 subs 和 vpa 常见的符号计算(六类运算) 因式分解、展开、合并、简化、通分和反函数等 计算极限 计算导数 计算积分 符号级数求和 代数方程和微分方程的求解(重点与难点)
因式分解 因式分解 factor(f) 例: factor 也可用于正整数的分解 例: syms x; f=x^6+1; factor(f) s=factor(100) 对大整数进行因式分解时可以先将其转化成符号常量 factor(12345678901234567890) % ERROR factor(sym('12345678901234567890'))
函数简化 函数简化 y=simplify(f) 例: 对符号表达式 f 进行简化 syms x; f=sin(x)^2 + cos(x)^2; y=simplify(f)
函数简化 函数简化 y=simple(f) 例:化简 对 f 尝试多种不同的方法(包括 simplify)进行简化, 以寻求其最简短形式 syms x; f=(cos(x)^2-sin(x)^2)*sin(2*x)*(exp(2*x) ... -2*exp(x)+1)/(exp(2*x)-1); y1=simplify(f) y2=simple(f)
计算极限 计算极限 limit(f,x,a) % 计算 limit(f,a) % 当默认变量趋向于 a 时的极限 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=limit((1-x/n)^n,n,inf)
计算导数 计算导数 g=diff(f,v) % 求符号表达式 f 关于变量 v 的导数 g=diff(f) % 计算关于默认变量的导数 g=diff(f,v,n) % 求 f 关于 v 的 n 阶导数 例: syms x; f=sin(x)+3*x^2; g1=diff(f,x) g2=diff(f,x,3)
计算积分 计算积分 int(f,v,a,b) % 计算定积分 int(f,a,b) % 计算关于默认变量的定积分 例:计算 和 syms x; f=(x^2+1)/(x^2-2*x+2)^2; I=int(f,x) K=int(exp(-x^2),x,0,inf)
符号级数求和 符号级数求和 symsum(f,v,a,b) % 级数求和 symsum(f,a,b) % 关于默认变量求和 syms n; f=1/n^2; S=symsum(f,n,1,inf) S100=symsum(f,n,1,100) 例:计算函数级数 syms n x; f=x/n^2; S=symsum(f,n,1,inf)
代数方程求解 代数方程求解 solve(f,v) % 求方程关于指定自变量的解 例: 这里 f 可以用字符串表示或符号表达式 得不到解析解时,给出数值解 例: solve('2*x-3') % 或 solve('2*x-3=0') syms x; solve(2*x-3) % 不能写成 solve(2*x-3=0) syms x; solve(2*x-sin(x)+1)
微分方程求解 ode23s、ode23tb 用 Maltab自带函数 解初值问题 求微分方程解析解:dsolve 求微分方程数值解*(自学,选学): ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb
dsolve 求解析解 dsolve 的使用 例 1:求微分方程 的通解,并验证。 y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v') 其中 y 为输出, eq1、eq2、...为微分方程,cond1、cond2、...为初值条件,v 为自变量。 例 1:求微分方程 的通解,并验证。 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') syms x; diff(y)+2*x*y - x*exp(-x^2)
dsolve 的使用 几点说明 微分方程中用 D 表示对 自变量 的导数,如: Dy y'; D2y y''; D3y y''' 如果省略初值条件,则表示求通解; 如果省略自变量,则默认自变量为 t dsolve('Dy=2*x','x'); % dy/dx = 2x dsolve('Dy=2*x'); % dy/dt = 2x 若找不到解析解,则返回其积分形式。
dsolve 的使用 例 2:求微分方程 在初值条件 下的特解,并画出解函数的图形。 例 2:求微分方程 在初值条件 下的特解,并画出解函数的图形。 y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y);
dsolve 的使用 例3:求微分方程组 在初值条件 下的特解,并画出解函数的图形。 例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]); 注:解微分方程组时,如果所给的输出个数与方程个数相同,则方程组的解按词典顺序输出;如果只给一个输出,则输出的是一个包含解的结构(structure)类型的数据。
dsolve 的使用 例: [x,y]=dsolve('Dx+5*x=0','Dy-3*y=0', ... 'x(0)=1', 'y(0)=1','t') r = dsolve('Dx+5*x=0','Dy-3*y=0', ... 'x(0)=1', 'y(0)=1','t') 这里返回的 r 是一个 结构类型 的数据 r.x %查看解函数 x(t) r.y %查看解函数 y(t) dsolve的输出个数只能为一个 或 与方程个数相等 只有很少一部分微分方程(组)能求出解析解。 大部分微分方程(组)只能利用数值方法求数值解。
求微分方程数值解*(自学,选学) [T,Y] = solver(odefun,tspan,y0) 其中 y0 为初值条件,tspan为求解区间;Matlab在数值求解时自动对求解区间进行分割,T (列向量) 中返回的是分割点的值(自变量),Y (数组) 中返回的是这些分割点上的近似解,其列数等于因变量的个数。 solver 为Matlab的ODE求解器(可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb) 没有一种算法可以有效地解决所有的 ODE 问题,因此MATLAB 提供了多种ODE求解器,对于不同的ODE,可以调用不同的求解器。
Matlab的ODE求解器*(自学,选学) 特点 说明 ode45 非刚性 单步法;4,5 阶 R-K 方法;累计截断误差为 (△x)3 大部分场合的首选方法 ode23 单步法;2,3 阶 R-K 方法;累计截断误差为 (△x)3 使用于精度较低的情形 ode113 多步法;Adams算法;高低精度均可到 10-3~10-6 计算时间比 ode45 短 ode23t 适度刚性 采用梯形算法 适度刚性情形 ode15s 刚性 多步法;Gear’s 反向数值微分;精度中等 若 ode45 失效时,可尝试使用 ode23s 单步法;2 阶Rosebrock 算法;低精度 当精度较低时,计算时间比 ode15s 短 ode23tb 梯形算法;低精度 当精度较低时,计算时间比ode15s短
参数说明*(自学,选学) [T,Y] = solver(odefun,tspan,y0) 例 : odefun 为显式常微分方程,可以用命令 inline 定义,或在函数文件中定义,然后通过函数句柄调用。 求初值问题 的数值解,求解范围为 [0,0.5] 例 : fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); 注:也可以在 tspan 中指定对求解区间的分割,如: [x,y]=ode23(fun,[0:0.1:0.5],1); %此时 x=[0:0.1:0.5]
上机作业 教材P28的习题1, 2, 3, 4 Matlab 演示