彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
4/1
线性方程组的解析解
符号积分变换
求解一元二次方程a * x 2 + b * x + c = 0的解。 ¿¿ f=sym(’a*x^2+b*x+c=0’); ¿¿ xf=solve(f)
彭朝霞 (北京航空航天大学)
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
11 / 1
微分方程解析解
微分方程的解析解dsolve输出顺序
求微分方程组 du =v dt dv =w dt dw = −u dt
并且初值满足： u (0) = 0, v (0) = 0, w (0) = 1. 通解的MATLAB程序为： ¿¿S1 = dsolve(’Du=v, Dv=w, Dw=-u’,’u(0)=0, v(0)=0, w(0)=1’) ¿¿S2 = dsolve(’Dw=-u,Du=v, Dv=w’, ’w(0)=1, u(0)=0, v(0)=0’) ¿¿S3 = dsolve(’Dw=-u,Dv=w, Du=v’, ’w(0)=1, v(0)=0, u(0)=0’) 注： 解微分方程组时， 如果所给的输出个数与方程个数相同， 则方程组的解按词典 顺序输出；
彭朝霞 (北京航空航天大学)
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
6/1
线性方程组的解析解
求下列线性代数方程组的解。
x1 cos(sita) − x2 sin(sita) = a x1 sin(sita) + x2 cos(sita) = b
¿¿ syms x1 x2 a b sita; ¿¿ L1=x1 * cos(sita)-x2 * sin(sita)-a; ¿¿ L2=x1 * sin(sita)+x2 * cos(sita)-b; ¿¿ [x1,x2]=solve(L1,L2,x1,x2) 运行结果： x1 = cos(sita)*a+b*sin(sita) x2 = -a*sin(sita)+cos(sita)*b
1
微分方程中用D 表示对自变量的导数/微分， 如： 求高阶微分： Dy → y ′ ; D 2y → y ′′ ; D 3y → y ′′′ 。任何D 后面所跟的字母为因变量。
2 3
如果省略初值条件， 则表示求通解； 如果省略自变量， 则默认自变量为t dsolve(’Dy=2*x’,’x’); % dy/dx = 2x dsolve(’Dy=2*x’); % dy/dt = 2x
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
8/1Biblioteka 微分方程解析解微分方程的解析解
函数dsolve用来解符号微分方程、 方程组。 dsolve 的使用 y=dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
如果没有初始条件， 则求出通解， 如果有初始条件， 则求出特解。
14 / 1
微分方程数值解
数值解法的基本思想
xi +1 − xi = h, i = 0, 1, 2, . . . , n − 1,用离散的方法求解微分方程： y ˙ = f (x , y ) y (x0 ) = y0 ♣ 欧拉方法,即用差商代替导数 ♣ 改进的欧拉方法 ♣ 使用泰勒公式 ♣ 泰勒公式法： 龙格-塔库法、 线性多步法等方法
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
3/1
线性方程组的解析解
代数方程的求解solve
解方程x 3 − 3 * x + 1 = 0. f 是符号表达式 ¿¿ syms x; f=x^3-3*x+1; s=solve(f,x) 注意： 不能写成s=solve(f=0,x) f 是字符串 ¿¿ s=solve(’x^3-3*x+1’,’x’) ¿¿ s=solve(’x^3-3*x+1=0’,’x’) 注意： 可以不含等号， 表示解方程f=0， 但是也可以含等号。
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
13 / 1
微分方程数值解
2、 微分方程的数值解
♣ 常微分方程数值解的定义: 在生产和科研中所处理的微分方程往往很复杂且大多得不出一般解.而在实际上对 初值问题， 一般是要求得到解在若干个点上满足规定精确度的近似值， 或者得到一 个满足精确度要求的便于计算的表达式. 因此， 研究常微分方程的数值解法是十分必要的. 对常微分方程 y ˙ = f (x , y ), y (x0 ) = y0
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
17 / 1
微分方程数值解
微分方程求数值解
第二步： 解微分方程 tfinal=100; x0=[0 0 1e-10]; [T , X ]=ode23(’zmfun’,[0,tfinal],x0)； plot(t,x)
MATLAB基础及工程应用
2014 年 10 月 23 日
5/1
线性方程组的解析解
代数方程求解
solve 也可以用来解方程组 solve (’eqn1’,’eqn2’,…,’eqnN’,’var1,var2,…,varN’) 求解由’eqn1’,’eqn2’,…,’eqnN’ 组成的方程组关于变量’var1,var2,…,varN’ 的解例 如： ¿¿ [x,y,z]=solve(’x+2*y-z=27’,’x+z=3’, ... ’x^2+3*y^2=28’,’x’,’y’,’z’) 输出变量的顺序要书写正确！ solve 在得不到解析解时， 会给出数值解
MATLAB基础及工程应用
2014 年 10 月 23 日
2/1
线性方程组的解析解
线性方程组的解析解
♣ 代数方程（非线性）求解： solve
s=solve(f,v) ： 求方程关于指定自变量的解 s=solve(f) ： 求方程关于默认自变量的解
f 可以是用字符串， 或符号表达式 来表示的方程
若f 是字符串， 可以不含等号， 表示解方程f=0 若f 是符号表达式， 不能含等号
[X,Y]=dsolve(’Dx+5*x+y=exp(t)’,’Dy-x-3*y=0’, ... ’x(0)=1’, ’y(0)=0’, ’t’) X =exp(2*t)/6 + (2*exp(t))/11 - 4*C3*exp(15^(1/2)*t - t) (4*C4)/exp(t + 15^(1/2)*t) + 15^(1/2)*C3*exp(15^(1/2)*t - t) (15^(1/2)*C4)/exp(t + 15^(1/2)*t) Y =C3*exp(15^(1/2)*t - t) - exp(t)/11 - (7*exp(2*t))/6 + C4/exp(t + 15^(1/2)*t) 注： 解微分方程组时， 如果所给的输出个数与方程个数相同， 则方程组的解按词典 顺序输出； 如果只给一个输出， 则输出的是一个包含解的结构（structure）类型的 数据。
第四章MATLAB 符号运算
彭朝霞 Email: pengzhaoxia@ 北京航空航天大学 交通科学与工程学院 2014 年 10 月 23 日
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
1/1
本章主要内容
彭朝霞 (北京航空航天大学)
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
15 / 1
微分方程数值解
2、 微分方程的数值解
♣ 微分方程数值求解函数与基本应用格式为： [t,y]=ode23(‘fun’,T0,Tfinal,Y0) [t,y]=ode45(‘fun’,T0,Tfinal,Y0) 函数说明： (1) ode23采用龙格-库塔二阶算法 ode45采用龙格-库塔四阶算法 (2) fun为由m函数定义的微分方程组的m函数名。 (3) T0,Tfinal,为数值积分区间， 仿真范围 (4) Y0为微分方程的初始条件向量。
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
7/1
线性方程组的解析解
求下列非线性代数方程组以y、 z 作变量的解。
uy 2 + vz + w = 0 y +z +w =0 ¿¿ syms y z u v w; ¿¿ eq1=u*y^2+v*z+w; ¿¿ eq2=y+z+w; ¿¿ [y z]=solve(eq1,eq2, y, z) y = [-1/2/u*(-2*u*w-v+(4*u*w*v+v^2-4*u*w)^(1/2))-w] [-1/2/u*(-2*u*w-v-(4*u*w*v+v^2-4*u*w)^(1/2))-w] z = [1/2/u*(-2*u*w-v+(4*u*w*v+v^2-4*u*w)^(1/2))] [1/2/u*(-2*u*w-v-(4*u*w*v+v^2-4*u*w)^(1/2))]
彭朝霞 (北京航空航天大学)
MATLAB基础及工程应用
2014 年 10 月 23 日
16 / 1
微分方程数值解
微分方程求数值解
例： 解下列微分方程组： 8 x ˙ 1 = − x1 + x2 x3 3 x ˙ 2 = −10x2 + 10x3 x ˙ 3 = −x1 x2 + 28x2 − x3 初始条件： x1 (0) = x2 (0) = 0, x3 = ������ = 1−10 。 1、 求微分方程的数值解 第一步： 建立M函数zmfum.m function dx=zmfun(t,x) dx=zeros(3,1); dx(1)=-8*x(1)/3+x(2)*x(3); dx(2)= -10*x(2)+10*x(3); dx(3)= -x(1)*x(2)+28*x(2)-x(3);