代数方程和微分方程求解.ppt

第4章 代数方程和微分方程求解 代数方程组的符号解 代数方程组的符号求解函数为solve函数 函数调用格式为x solve fun1 fun2 funm var1 varm 其中fun1 fun2 为符号表达式 var1 var2 为符号变量 指明所求解的变量 如果不指明所求解的变量 则系统根据变量在英文字母中与x的接近程度选择 例4 1求解线性方程组 参考程序symsxyza u v w solve a x y z 1 x a y z 2 x y a z 3 x y z 输出结果 u a 4 a 2 a 2 v 2 a 2 w 3 a a 2 a 2 也可以利用下面的语句求解 s solve a x y z 1 x a y z 2 x y a z 3 x y z s x 1x1sym y 1x1sym z 1x1sym 其中 s是结构形数据 可以用下面的方法显示求得的解 s xans a 4 a 2 a 2 方程组的符号求解得到的是解析式 其优点是可以得到含未知参数的解 有利于进一步的分析 缺点是计算效率低 且复杂的问题一般得不到解析解 常微分方程组的符号求解 常微分方程 组 的符号运算函数为dsolve 调用格式为变量 dsolve eq1 eq2 其中equ1 equ2等是字符串 描述微分方程 方程形式与solve语句中类似 其中的导数项利用D表示 Dy D2y分别描述y的一阶和二阶导数 初值和边值条件也利用方程描述 例4 2 求微分方程的通解和在初始条件y 0 1 y 0 2下的特解 y dsolve D2y 2 Dy 2 y 0 y C24 exp t cos t C25 exp t sin t y dsolve D2y 2 Dy 2 y 0 y 0 1 Dy 0 2 y exp t cos t exp t sin t 代数方程的数值方法 多项式在matlab语言中 多项式利用行向量表示 如向量u 120 54 用作多项式函数时 与多项式对应 因此 多项式运算对应向量的相关运算 如多项式的加法对应向量的加法等 多项式求根 多项式运算的几个常用函数 P conv p1 p2 多项式乘法 d r deconv p1 p2 多项式除法Dp polyder p 多项式的导数Ip polyint p 多项式的积分 原函数 Y polyval p x 输出多项式p在向量x的值 多项式求根的函数为r roots p 求得多项式p的所有根 例4 3 求多项式的根并在多项式图形中表示 参考程序 q 1 20138 328 2231692 1260 r roots q x 2 2 0 05 8 y polyval q x y1 polyval q r plot x y r y1 p xlim 2 2 8 legend polynomial roots 线性方程组的求解 线性方程组Ax b可以利用矩阵除法直接得到 但当系数矩阵为稀疏矩阵时 利用稀疏矩阵函数可以得到更高的计算效率 稀疏矩阵利用函数A1 sparse A 定义 例4 4 求解n阶线性方程组对n 1000 分别利用正常方法和稀疏矩阵求解 并比较计算时间 参考程序 n 1000 A diag 2 ones 1 n diag ones 1 n 1 1 diag ones 1 n 1 1 A n 1 1 b ones n 1 ticx A b tocA1 sparse A ticx A1 b toc输出结果 Elapsedtimeis0 144504seconds Elapsedtimeis0 009240seconds 求解非线性方程f x 0一般利用以下形式的迭代方法 xk 1 xk k 1 2 如著名的Newton迭代法就是众多迭代方法中的一种 非线性方程求根 迭代是一种逐步近似的过程 在迭代求根中 需要给定初始近似来启动迭代过程 初始近似的选取非常重要 这是因为许多方法只有当接近所求的解时 迭代才收敛 当方程有多个解时 需要充分靠近所求的解 才能保证迭代收敛到这个解 Matlab中非线性方程求根的函数是fzero 调用格式为p fzero fun x0 其中x0是初始值 fun可以是字符串以指示函数名 也可以是以 开头的函数句柄如 sin sin x sin x 都可以 应用最灵活的是后一种 例4 5 求方程sin x cos x 1的解 初始值依次选取x0 1 2 3 10 参考程序 r zeros 1 10 forx0 1 10r x0 fzero x sin x cos x 1 x0 endr运行结果 r 1 57081 57081 57081 57086 28326 28326 28327 85407 85407 8540运行结果显示 在一般情况下 迭代收敛到与初始值接近的根 例4 6 求解方程 functionr rootexample x0 r fzero rootfun x0 endfunctionf rootfun x f x quad t exp t 2 0 x 1 end运行结果 r rootexample 0 8 r 0 5219 例4 7 求方程的最小正根 首先寻找初始值 画出方程在 0 2 上的图形 x 0 pi 50 2 pi y x sin x 2 x sin 2 x 3 x sin 3 x 4 x sin 4 x 5 x sin 5 x 10 plot x y x 0 利用对上述图形的局部放大图 下图 估计得到初始值3 8 程序和结果 formatlong r fzero x x sin x 2 x sin 2 x 3 x sin 3 x 4 x sin 4 x 5 x sin 5 x 10 3 8 r 3 908485977246559 例4 8 求含参数代数方程的根与参数a的关系 注 含参数的方程的求解对后续计算非常重要 如上面例子求参数 使得根取极大等 functionrootfindingx 0 0 01 5 l 1 root 1 zeros 1 length x fora x root 1 l fzero x fun x a a l l 1 endplot x root 1 xlabel parameters ylabel roots endfunctionf fun x a f a 3 x 2 exp 0 2 x exp 0 02 x sin a x 4 end 非线性方程组求解 非线性方程组的一般形式为F x 0其中 非线性方程组求解的matlab函数为x fsolve fun x0 使用方法完全类似于函数fzero 但对初始条件的选择要困难的多 例4 9 求解非线性方程组的全部实根 解 首先将两条联立的曲线画出 以确定根 交点 的个数和位置 参考程序 h ezplot x y x 2 cos pi x y 2 sin pi y pi 2 3 3 3 3 str1 get get gca Title String set h Color 100 LineStyle holdonh1 ezplot x y x 2 y 2 2 sin 2 x 2 y 2 4 3 3 3 3 str2 get get gca Title String title str1 str2 gridon 输出结果 由上图知 方程组共有20个不同的解 再由方程组关于x的对称性知只要求出右边的10个 利用局部放大图 如下图 可得到10个解的近似值 数值程序 functionr d2solveh ezplot x y x 2 cos pi x y 2 sin pi y pi 2 3 3 3 3 str1 get get gca Title String set h Color 100 LineStyle holdonh1 ezplot x y x 2 y 2 2 sin 2 x 2 y 2 4 3 3 3 3 str2 get get gca Title String title str1 str2 x 1 72 1 7 1 68 1 68 1 66 1 64 0 36 0 65 0 9 0 88 y 0 63 0 42 0 19 0 16 0 42 0 7 1 85 1 38 1 7 1 53 X x y m length x r zeros 2 m fork 1 mr k fsolve solfun X k endplot r 1 r 2 ko endfunctionf solfun x f x 1 2 cos pi x 1 x 2 2 sin pi x 2 pi 2x 1 2 x 2 2 2 sin 2 x 1 2 x 2 2 4 end 程序运行结果 r d2solver 1 72511 70381 68741 68521 66791 64300 36050 64870 89800 8899 0 6279 0 4121 0 14940 15040 41970 7085 1 8548 1 3778 1 7007 1 5239 常微分方程 组 的数值解 一阶常微分方程组的显式形式为 1 其中一般的高阶常微分方程组都可以通过将低阶导数作为新的变量化为上述形式 例4 10 考虑炮弹的抛射轨迹 在不考虑空气阻力时 利用Newton第二定律 其轨迹可以由以下微分方程组描述 取则方程组 4 2 化为其中 Matlab函数针对微分方程的一般形式 1 求解 初值问题 求解初值问题的matlab函数采用自适应Runge Kutta方法 即自变量节点的间距按精度调整 因此 自变量节点是不等距的 其中 fun是微分方程组 4 3 右端函数F 是自变量取值范围 x0是初值向量 t是输出的自变量节点向量 x的各列是未知函数各分量对应节点向量的数值解 较复杂的用法是这里 options是算法的有关参数的设置变量 它利用函数odeset赋值 a b 是方程中的参数 常用的两个函数是x ode45 fun t0 t1 x0 4 5阶Runge Kutta算法x ode23 fun t0 t1 x0 2 3阶Runge Kutta算法 odeset函数的调用格式为options odeset 属性名 属性值 例4 9 考虑炮弹的抛射轨迹方程 x 0 0 y 0 h x 0 v0cos y 0 v0sin 计算不同角度的投射轨道 functionproj thisfunctioncomputestheprojectsofabulletinspace Eachcurve describesamovingtrailwithdifferentshootingangle Theinitialheightoftheabulletish 12 Theinitialvelocityofthebulletisv0 20 h1 plot 0 60 0 0 b LineWidth 3 holdonh 12 v0 20 foralpha 0 pi 24 pi 2projmain h v0 alpha endend functionprojmain h v0 alpha X0 0 v0 cos alpha h v0 sin alpha options odeset Reltol 1e 10 t X od