数字信号处理课程设计
题目：基于MATLAB的线性常系数差分方程求解学院：
专业：
班级：
学号：
姓名：
指导教师：
目录
摘要 (1)
第一章背景 (3)
1.1 背景知识 (3)
1.2 《数字信号课程》特点 (3)
1.3 软件介绍 (4)
1.4 MATLAB及数字信号处理 (4)
第二章设计目的及要求 (6)
2.1 设计目的 (6)
2.2 课程设计的内容要求 (7)
2.2.1 设计要求 (7)
第三章设计任务 (8)
第四章设计原理 (9)
4.1 差分与差分方程 (9)
4.2 线性常系数差分方程 (14)
4.3 线性常系数差分方程的求解 (15)
第五章设计过程 (16)
5.1 用MATLAB求解差分方程 (16)
第六章设计代码及结果 (18)
6.1 MATLAB源程序 (18)
6.2 程序运行结果 (20)
6.3 比较结果总结 (24)
第七章收获与体会 (25)
致谢 (27)
参考文献 (28)
摘要
《数字信号处理》分析了数字信号处理课程的重要性及特点，为了帮助学生理解与掌握课程中的基本概念、基本原理、基本分析方法，提出了用MATLAB进行数字信号处理课程设计的思路，并阐述了课程设计的具体方法、步骤和内容。
MATLAB语言是一种广泛应用于工程计算及数值分析领域的新型高级语言，MATLAB功能强大、简单易学、编成效率高，深受广大科技工作者的喜爱，特别是MATLAB还具有信号分析工具箱，不需具备很强的编程能力，就可以很方便地进行语音信号分析、处理和设计。
线性常系数差分方程求解是数字信号处理课程中常出现的课题，也是现代科学中值得深入研究的一个课题
本文介绍了线性常系数差分方程的基本概念，论述了其求解方法，并用MATLAB具体实现了线性常系数差分方程的求解。
基于MATLAB的线性常系数差分方程求解主要是用MATLAB作为工具平台，设计中涉及到差分方程的递推求解以及用filter对系数向量的归一化等等。通过数字信号处理课程的理论知识的综合运用，从实践上初步实现对数字信号的处理。
关键字：MATLAB，线性常系数差分方程，数字信号处理。
Abstract
" Digital signal processing " analysis of digital signal processing course of the importance and features, in order to help the students to understand and grasp basic concepts, basic principles, basic analysis method, is put forward with the MATLAB for digital signal processing curriculum design, curriculum design and describes the specific methods, steps and content.
The MATLAB language is widely used in engineering calculation and the numerical analysis in the field of advanced language, MATLAB powerful, easy to learn, a high efficiency, by the vast number of scientific workers' favorite, especially the MATLAB also has a signal analysis toolbox, does not need to have very strong ability of programming, can be very convenient for the analysis of speech signal, processing and design.
Linear constant coefficient differential equation is a digital signal processing program that often appear in the topic, as well as modern science and worthy of in-depth study of a topic
This paper introduces the linear constant coefficient differential equation basic concept, discussed the solution method, and MATLAB the specific realization of linear constant coefficient difference equation.
MATLAB based on the linear constant coefficient difference equation MATLAB is mainly used as a tool platform, design relate to differential equation recursive solution and the use of filter on the coefficient vector is normalized and so on. Through the course of digital signal processing theory, the integrated use of knowledge, from the practice of preliminary implementation of digital signal processing.
Key words: MATLAB,Linear constant coefficient differential equation,Digital signal processing.
第一章背景
1.1 背景知识
数字信号处理（Digital Signal Processing,简称DSP）是一门设计许多学科而又广泛应用于许多领域的新兴学科。DSP有两种含义：Digital Signal Processing（数字信号处理）、Digital Signal Processor（数字信号处理器）。我们常说的
DSP指的是数字信号处理器。数字信号处理器是一种适合完成数字信号处理运算的处理器。20世纪60年代以来，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。在过去的二十多年的时间里，数字信号处理已经在通信等领域得到极为广泛的应用。
数字信号处理是利用计算机专用处理设备，以数字形式对信号进行采集、变换、滤波、估值增强、压缩、识别等处理，以得到符合人们需要的信号形式。它是以众多学科为理论基础的，它所涉及的范围极其广泛。例如，在数学领域，微积分、概率统计、随即过程、数值分析等都是数字信号处理的基本工具，与网络理论、信号与系统、控制论、通信理论、故障诊断等也密切相关，近来新兴的一些学科，如人工智能、模式识别、神经网络等，都与数字信号处理密不可分。可以说，数字信号处理是把许多经典的理论体系作为自己的理论基础，同时又使自己成为一系列新兴学科的理论基础。
1.2 《数字信号课程》特点
《数字信号处理》课程是一门理论和技术发展十分迅速、应用非常广泛的前沿性学科，他的理论性和实践性都很强，他的特点是：
(1)要求的数学知识多，包括高等代数、数值分析、概率统计、随机过程等。
(2)要求掌握的基础知识强，网络理论、信号与系统是本课程的理论基础。
(3)与其他学科密切相关，即与通信理论、计算机、微电子技术不可分，又是人工智能、模式识别、神经网络等新兴学科的理论基础之一。
选择用MATLAB进行课程设计:MATLAB语言是一种广泛应用于工程计算及数值分析领域的新型高级语言，MATLAB功能强大、简单易学、编程效率高，深受广大科技工作者的欢迎。特别是MATLAB还具有信号分析工具箱，不需具备很强的编程能力，就可以很方便地进行信号分析、处理和设计。
1.3 软件介绍
MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数字数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、经融建模设计与分析等领域。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解答问题要比用C,FORTRAN等语言完成相同的事情简捷的多，并且mathworks也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C,FORTRAN，C++，JAVA的支持。可以直接调用，用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。
1.4 MATLAB及数字信号处理
MATLAB是矩阵实验室之意。除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学，工程中常用的形式十分相似，故用MATLAB来解决问题要比用C,FORTRAN等语言完全相同的事情简捷得多，可以将自己编写的实用程序导入到MATLAB函数库中方便在新的版本中也加入了对C，FORTRAN，c++，JAVA的支持。可以直接调用，用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用。
信号时数字信号处理领域中最基本、最重要的概念。简单地说，信号就是信息的载体，是信息的物理体现。信号既可以分为时间连续、幅度也连续的模拟信号和时间和幅度上都经过量化的数字信号，也可以划分为连续时间信号和离散时间信号。几乎在科学技术的每一领域，为了信号的提取，都要进行信号处理，就是以数值计算的方法对信号进行采集，变换，综合，估计，与识别的加工处理过程，借以达到提取信息和便与应用的目的。随着计算机和信息科学的飞速发展，信号处理已经逐渐发展为一门独立的学科，是信息科学的重要组成部分。在语音处理，雷达，航空，航天，地质勘探，通信，生物医学工程等众多领域得到了广泛应用。
MATLAB软件，在数字信号处理方面具有得天独厚的优势。利用目录下的系统函数，用户可以实现波形的产生，信号的变换，滤波，功率谱估计，系统设计
与稳定性分析，小波信号分析等众多功能。本文即是以数字信号处理的理论基础，应用MATLAB软件求解线性常系数差分方程的一个具体事例。
第二章设计目的及要求
2.1 设计目的
1.全面复习课程所学理论知识，巩固所学知识重点和难点，将理论与实践
很好的结合起来。
2.掌握信号分析与处理的基本方法与实现
3.提高综合运用所学知识独立分析和解决问题的能力
4.熟练使用一种高级语言进行编程实现
课程设计是教学的最后一个步骤，课程设计有利于基础知识的理解，我们掌握了基础知识和基本技能，但是要真正接触才能真正理解课程的深入部分；还有礼物逻辑思维的锻炼，在许多常规学科的日常教学中，我们不难发现这样一个现象，不少学生的思维常常处于混乱的状态，写起作文来前言不搭后语，解起数学题来步骤混乱，这些都是缺乏思维训练的结果，所以我们可以通过实践来分析问题、解决问题、预测目标等；同时也有利于其他学科的整合，例如我们这次的课程设计就是要运用MATLAB软件的帮助才能实现；最重要的有利于治学态度的培养，在课程设计中，我们可能经常犯很多小错误，可能要通过好几次的反复修改、调试才能成功，但这种现象会随着学习的深入而慢慢改观。这当中就有一个严谨治学、一丝不苟的科学精神的培养，又有一个不拍失败、百折不挠品格的锻炼。
《数字信号处理》课程设计实在学生完成数字信号处理和MATLAB的结合后基于实验以后开设的。本课程设计的目的是为了让学生综合数字信号处理和MATLAB并实现线性常系数差分方程求解。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序编写与调试、写课程设计报告等步骤，使学生进一步掌握数字信号处理课程的基本理论、基本方法和基本技术；使学生增进对MATLAB 的认识，利用MATLAB加深对理论知识的理解；使学生了解和掌握使用MATLAB
的应用过程和方法，为以后的设计打下良好基础培养学生能根据设计要求，进行理论知识分析、设计方法总结、典型实例设计等方面的设计综合能力；使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。
随着信息科学和计算机技术的迅速发展，现代信号处理的理论与应用得到飞跃式的发展，形成一门极其重要的学科。特别是数字信号处理也已经成为高等学校相关专业的必修课程。为了更好地将数字信号处理的理论付诸实践，此次课程设计是一个很好的契机，通过仿真实验，让大家初步了解信号处理的分析法和实现方法。
现代信号处理是一门一算法为核心，理论和实践性较强的学科。是电子信息工程、通信工程专业、电子信息科学与技术专业的一门重要的专业基础课。数字信号处理课程是在学习完数字信号处理的相关理论后，进行的综合性训练课程，其目的是：使学生进一步巩固数字信号处理的基本盖帘、理论、分析方法和实现
方法；增强学生应用MATLAB语言编写数字信号处理的应用程序及分析、解决实际问题的能力。
2.2 课程设计的内容要求
2.2.1 设计要求
1.巩固所学的专业技术知识；
2.提高综合运用所学理论知识独立分析和解决问题的能力；
3.进一步提高程序设计及调试能力；
4.更好地将理论与实践相结合；
5.学习和掌握科学研究资料检索的方法，学习对已有资料进行消化总结的方法
6.学习撰写科学报告的基本方法；
7.本设计要求分组合作完成；
8.上机前提前熟悉使用课程设计实验平台，掌握其基本的操作方法；
9.上机前了解课程设计实验平台的源代码，掌握其程序结构及在此平台上添加
处理程序的方法；
10.设计过程中详细记录产生的图形、参数、数据等，用于编写课程设计报告。
本次课程设计的部分课题用到了本科数字信号处理的几乎所有知识，既能帮助学生对知识点的融会贯通，又使学生感到学有所用，培养进一步深入学习信号处理的兴趣。
2.2.2 设计内容
1.课程设计题目和题目涉及要求
2.设计思想和系统分析功能分析
3.设计中关键部分的理论分析与计算，关键模块的设计思路
4.程序代码清单
5.测试数据、测试输出结果，及必要的理论分析和比较
6.总结，包括设计过程中遇到的问题和解决办法，设计心得与体会等
7.致谢
8.参考文献
题目：基于MATLAB 的线性常系数差分方程求解
1、自行产生一个序列，要求：
（1）对序列进行差分运算，并画出差分序列的时域波形图；
2、已知一个二阶线性常系数差分方程用下式表示：
y(n)+a 1y(n-1)+a 2y(n-2)= b 0x(n)+b 1x(n-1)+b 2x(n-2)，
要求：
（1） 参数a 1、a 2、b 0、b 1、b 2由运行时输入；
（2） 已知输入)(5.0)(n u n x n ，画出x(n)的时域波形图；
（3） 求出x(n)的共轭对称分量x e (n)和共轭反对称分量x o (n)，并分别画出时
域波形图；
（4） 初始条件由运行时输入，求输出y(n)，并画出其波形；
（5） 对于不同的初始条件，分析其输出是否一致，从中得出什么结论
4.1 差分与差分方程
与连续时间信号的微分及积分运算相对应，离散时间信号有差分及序列求和运算。设有序列f(k),则称…,f(k+2),f(k+1),…,f(k －1),f(k －2),…为f(k)的移位序列。序列的差分可以分为前向差分和后向差分。一阶前向差分定义为
()(1)()f k f k f k ?=+- （3.1—1）
一阶后向差分定义为
()()(1)f k f k f k ?=-- （3.1—2）
式中Δ和Δ称为差分算子。由式（3.1—1）和式（3.1—2）可见，前向差分与后向差分的关系为
()(1)f k f k ?=?- （3.1—3）
二者仅移位不同，没有原则上的差别，因而它们的性质也相同。此处主要采用后向差分，并简称其为差分。
由查分的定义，若有序列1()f k 、2()f k 和常数1a ，2a 则
1122112211221112221122[()()][()()][(1)(1)]
[()(1)][()(1)]
()()
a f k a f k a f k a f k a f k a f k a f k f k a f k f k a f k a f k ?+=+--+-=--+--=?+? （3.1—4）
这表明差分运算具有线性性质。
二阶差分可定义为 2()[()][()(1)]()(1)
()2(1)(2)f k f k f k f k f k f k f k f k f k ?=??=?--=?-?-=--+- （3.1—5）
类似的，可定义三阶、四阶、…、n 阶差分。一般地，n 阶差分
10()[()](1)()n
n n j j n f k f k f k j j -=???=??=-- ???∑ （3.1—6）
式中
!,0,1,2,,()!!
n n j n j n j j ??== ?-?? （3.1—7）
为二项式系数
序列f(k)的求和运算为 ()k
i f i =-∝∑ （3.1—8） 差分方程是包含关于变量k 的未知序列y(k)及其各阶差分的方程式，它的一般形式可写为
,(),(),,()0n F k y k y k y k ????=?? （3.1—9a ）
式中差分的最高阶为n 阶，称为n 阶差分方程。由式（3.1—6）可知，各阶差分均可写为y(k)及其各移位序列的线性组合，故上式常写为
[],(),(1),,()0G k y k y k y k n --= （3.1—9b ）
通常所说的差分方程是指式（3.1—9b ）形式的方程。
若式（3.1—9b ）中，y(k)及其各移位序列均为常数，就称其为常系数差分方程；如果某些系数是变量k 的函数，就称其为变系数差分方程。描述LTI 离散系统的是常系数线性差分方程。
差分方程是具有递推关系的代数方程，若一直初始条件和激励，利用迭代法渴求的差分方程的数值解。
4.1.1 差分方程的经典解
一般而言，如果但输入—单输出的LTI 系统的激励f(k)，其全响应为y(k)，那么，描述该系统激励f(k)与响应y(k)之间关系的数学模型式n 阶常系数线性差分方程，它可写为
1010()(1)()
()(1)()n m m y k a y k a y k n b f k b f k b f k m --+-++-=+-+
+- （3.1—10a ） 式中(0,1,,1)j a j n =-、(0,1,,)i b i m =都是常数。上式可缩写为
00
()()(=1)n m n j m i n j i a
y k j a f k i a --==-=-∑∑式中 （3.1—10b ） 与微分方程的经典解类似，上述差分方程的解由齐次解和特解两部分组成。齐次解用()h y k 表示，特解用()p y k 表示，即
()=()()h p y k y k y k + （3.1—11）
a.齐次解
当式（3.1—10)中的f(k)及其各移位项均为零时，齐次方程
10()(1)()0n y k a y k a y k n -+-+
+-= （3.1—12）
的解称为齐次解。 首先分析最简单的一阶差分方程。若一阶差分方程的齐次方程为
()(1)0y k ay k +-= （3.1—13）
它可改写为
()(1)y k a y k =--
y(k)与y(k －1)之比等于－a 表明，序列y(k)是一个公比为－a 的等比级数，因此y(k)应有如下形式
()()k
y k C a =- （3.1—14）
式中C 式常数，有初始条件确定。
对于n 阶齐次差分方程，它的齐次解由形式为k C λ的序列组合而成，将
k C λ代入到式（3.1—12），得
111100k k k n k n n C a C a C a C λλλλ-----++
++= 由于C ≠0，消去C ；且λ≠0，以k n λ-除上式，得
11100n n n a a a λλλ--++++=（3.1—15）
上式称为差分方程式（3.1—10）和式（3.1—12）的特征方程，它有n 个根
(0,1,,)j j n λ=，称为差分方程的特征根。显然，形式为j
k j C λ的序列都满足式
（3.1—12），因而它们是式（3.1—10）方程的齐次解。依特征根取值的不同，
差分方程齐次解的形式见表3—1,其中j C 、j D 、j A 、j θ等为待定常数
表3—1 不同特征根所对应的齐次解 sin()]k A βρ或0)]θ-b.特解
特解的函数形式与激励的函数形式有关，表3—2列出了集中典型的激励f(k)所对应的特解()p y k 。选定特解后代入原差分方程，求出其待定系数
()j P A θ或等，就得出方程的特解。
表3—2 不同激励所对应的特解 ),A j e θ其中
c.全解
式（3.1—10）的线性差分方程的全解是齐次解与特解之和。如果方程的特征根均为单根，则差分方程的全解为
1()=()()()
n k h p j j p j y k y k y k C y k λ=+=+∑（3.1—16）
如果特征根1λ为r 重根，而其余n －r 个特征根为单根时，差分方程的全解为
11()=()r n r j k k j j j j p j j r y k C k
C y k λλ-==+++∑∑（3.1—17）
式中各系数
j C 由初始条件确定。
如果激励信号是在k=0时接入的，差分方程的解适合于k ≥0。对于n 阶差分方程，用给定的n 个初始条件y(0),y(1),…,y(n －1)就可确定全部待定系数。如果差分方程的特解都是单根，则方程的全解为式（3.1—16），将给定的初始条件y(0),y(1),…,y(n －1)分别代入到式（3.1—16），可得
1211221111122(0)(0)
(1)(1)
(1)(1)n p n n p n n n n n p y C C C y y C C C y y n C C C y n λλλλλλ
---=++++=++++-=++
++-（3.1—18） 由以上方程可求得全部待定系数(0,1,,)j C j
n =。
4.1.2 零输入响应 系统的激励为零，仅由系统的初始状态引起的响应，称为零输入响应，用()zi y k 表示。在零输入条件下，式（3.1—10）等号右端为零，化为齐次方程，即
0()0n n j zi j a
y k j -=-=∑ （3.1—25）
一般设定激励是在k=0时接入系统的，在k ＜0时，激励尚未接入，故式（3.1—25）的几个初始状态满足
(1)(1)
(2)(2)
()()
zi zi zi y y y y y n y n -=--=--=
- （3.1—26）
式（3.1—26）中的y(－1),y(－2),…,y(－n)为系数的初始状态，由式（3.1—25）和式（3.1—26）可求得零输入响应
()zi y k 。
4.1.3 零状态响应
当系统的初始状态为零，仅由激励f(k)所产生的响应，称为零状态响应，用 表示。在零状态情况下，式（3.1—10）仍是非齐次方程，其初始状态为零，即零状态响应满足 00()()(1)(2)()0n m
n j zs m i j i zs zs zs a
y k j a f k i y y y n --==-=--=-==-=∑∑（3.1—30）
的解。若其特征根均为单根，则其零状态响应为
1()()
n k zs zsj j p j y k C y k λ==+∑ （3.1—31）
式中zsj C 为待定常数，()p y k 为特解。需要指出，零状态响应的初始状态
(1),(2),,()zs zs zs y y y n ---为零为零，但其初始值(0),(1),,(1)zs zs zs y y y n -不一定等于零。
4.2 线性常系数差分方程
4.2．1 一个N 阶线性常系数差分方程可用下式表示：
01()()()
M N
i i i i y n b x n i a y n i ===---∑∑
（1.4.1）
或者
000()(),1N M i i
i i a y n i b x n i a ==-=-=∑∑ （1.4.2）
式中，x （n ）和y(n)分别是系统的输入序列和输出序列，ai 和bi 均为常系数，式中y(n-i)和x(n-i)项只有一次幂，也没有相互交叉相乘项，故称为线性常系数差分方程。差分方程的阶数是用方程y(n-i)项中i 的最大取值与最小取值之差确定的。在（1.4.2）式中，y(n-i)项i 最大的取值N ，i 的最小取值为零，因此称为N 阶差分方程。
4.3 线性常系数差分方程的求解
已知系统的输入序列，通过求解差分方程可以求出输出序列。求解差分方程的基本方法有以下三种：
（1） 经典解法。这种方法类似于模拟系统中求解微分方程的方法，它包括齐次解与
特解，由边界条件求待定系数，上节已作简单介绍，这里不作介绍。
（2） 递推解法。这种方法简单，且适合用计算机求解，但只能得到数值解，对于阶
次较高的线性常系数差分方程不容易得到封闭式（公式）解答。
（3） 变换域方法。这种方法是将差分方程变换到z 域进行求解，方法简便有效。
当然还可以不直接求解差分方程，而是先由差分方程求出系统的单位脉冲响应，再与已知的输入序列进行卷积运算，得到系统输出。但是系统的单位脉冲响应如果不是预先知道，仍然需要求解差分方程，求其零状态响应解。
（4） 卷积法：由差分方程求出系统的h(n)，再与已知的x(n) 进行卷积，得到
y(n)。
观察（1.4.1）式，求n 时刻的输出，要知道n 时刻以及n 时刻以前的输入序列值，还要知道n 时刻以前的N 个输出序列值。因此求解差分方程在给定输入序列的条件下，还需要确定N 个初始条件。如果求n0时刻以后的输出，n0时刻以前N 个输出值y （n0-1）、y （n0-2）、??????、y （n0-N ）就构成了初始条件。 （1.4.1）式表明，已知输入序列和N 个初始条件，则可以求出n 时刻的输出；如果将该公式中的n 用n+1代替，可以求出n+1时刻的输出，因此（1.4.1） 式表示的差分方程本身就是一个适合递推法求解的方程。
第五章设计过程
5.1 用MATLAB求解差分方程
MATLAB信号处理工具箱提供的filter函数实现线性常系数差分方程的递推求解，调用格式如下：
yn=filter(B,A.xn) 计算系统对输入信号向量xn的零状态响应输出信号向量yn，yn与xn长度相等，其中，B和A是（1.4.2）式所给差分方程的系数向量，即
B=[b0,b1,???,bM], A=[a0,a1,???,aN]
其中a0=1，如果a0≠1，则filter用a0对系数向量B和A归一化。
yn=filter(B,A.xn,xi) 计算系统对输入信号向量xn的全响应输出信号yn。所谓全响应，就是由初始状态引起的零输入响应和由输入信号xn引起的零状态响应之和。其中，xi是等效初始条件的输入序列，所以xi是由初始条件确定的。MATLAB信号处理工具箱提供的filtic就是由初始条件计算xi的函数，其调用格式如下：
xi=filtic（B,A,ys,xs）
其中，ys和xs是初始条件向量：ys=[y(-1),y(-2),y(-3),??
?,y(-N)],xs=[x(-1),x(-2),x(-3),???,x(-M)]。如果xn是因果序列，则xs=0.调用时可缺省xs。
例1.4.1的MATLAB求解程序ep141.m如下：
%1.4.1.m:调用MATLAB解差分方程y(n)-0.8y(n-1)=x(n)
a=0.8；ys=1; %设差分方程系数a=0.8，初始状态：y(-1)=1
xn=[1,zeros(1,30)]; %x(n)=单位脉冲序列，长度N=31
B=1;A=[1,-0.8]; %差分方程系数
xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xi
yn=filter(B,A,xn,xi); %调用fiter解差分方程，求系统输出信号y(n)
n=0:length(yn)-1;
stem(n,yn,'.')
title('时域波形图)');xlabel('n');ylabel('y(n)')
程序中取查分方程系数a=0.8时，得到系统输出y（n）如图1.4.1（a）所示，与例1.4.1的解析递推结果完全相同。如果令初始条件y（-1）=0（仅修改程序中ys=0），则得到系统输出y（n）=h（n），如图1.4.1（b）所示。
(a)(b)
图（a）为a=0.8，y（-1）=1时，系统输出时域波形图，图（b）为a=0.8，y(-1)=0时，系统输出时域波形图。
第六章设计代码及结果
6.1 MATLAB源程序
6.1.1 源程序如下：
%1.m:调用MATLAB解差分方程y(n)-0.8y(n-1)=x(n)
ys=1; %初始状态：y(-1)=1
xn=[1,zeros(1,30)]; %x(n)=单位脉冲序列，长度N=31
B=1;A=[1,-0.8]; %差分方程系数
xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xi
yn=filter(B,A,xn,xi); %调用filter解差分方程，求系统输出信号y(n)
n=0:length(yn)-1; %n的取值范围
stem(n,yn,'.') %画出时域波形图
title('时域波形图)');xlabel('n');ylabel('y(n)') %x轴、y轴分别代表n，x(n)
6.1.2-2 源程序如下：
n=-5:5; %n的取值范围
xn=0.5.^n; %xn=0.5.^n
stem(n,xn,'fill'),grid on %画出时域波形图
xlabel('n'),ylabel('x(n)'),title('时域波形图') %x轴、y轴分别代表n，x(n)
6.1.2-3 源程序如下：
n=-5:5; %n的取值范围
a=0.5; %设a=0.5
Broyden方法求解非线性方程组的Matlab实现 注：matlab代码来自网络，仅供学习参考。 1.把以下代码复制在一个.m文件上 function [sol, it_hist, ierr] = brsola(x,f,tol, parms) % Broyden's Method solver, globally convergent % solver for f(x) = 0, Armijo rule, one vector storage % % This code es with no guarantee or warranty of any kind. % % function [sol, it_hist, ierr] = brsola(x,f,tol,parms) % % inputs: % initial iterate = x % function = f % tol = [atol, rtol] relative/absolute % error tolerances for the nonlinear iteration % parms = [maxit, maxdim] % maxit = maxmium number of nonlinear iterations % default = 40 % maxdim = maximum number of Broyden iterations % before restart, so maxdim-1 vectors are % stored % default = 40 % % output: % sol = solution % it_hist(maxit,3) = scaled l2 norms of nonlinear residuals % for the iteration, number function evaluations, % and number of steplength reductions % ierr = 0 upon successful termination % ierr = 1 if after maxit iterations % the termination criterion is not satsified. % ierr = 2 failure in the line search. The iteration % is terminated if too many steplength reductions % are taken. % % % internal parameter: % debug = turns on/off iteration statistics display as % the iteration progresses
解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;
if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵
MatLab解线性方程组一文通 当齐次线性方程AX=0,rank(A)=r 在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法. 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB 求解线性方程组 solve，linsolve 例： A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1]; %矩阵的行之间用分号隔开，元素之间用逗号或空格 B=[3;1;1;0] X=zeros(4,1);%建立一个4元列向量 X=linsolve(A,B) diff（fun，var，n）：对表达式fun中的变量var求n阶导数。 例如：F=sym（'u(x,y)*v(x,y)'）; %sym（）用来定义一个符号表达式 diff(F); %matlab区分大小写 pretty(ans) %pretty（）：用习惯书写方式显示变量；ans是答案表达式 非线性方程求解 fsolve(fun,x0,options) 为待解方程或方程组的文件名；fun其中 x0位求解方程的初始向量或矩阵； option为设置命令参数 建立文件fun.m： function y=fun(x) y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ... x(2) - 0.5*cos(x(1))+0.3*sin(x(2))]; >>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve')) 注： ...为续行符 m文件必须以function为文件头，调用符为@；文件名必须与定义的函数名相同；fsolve（）主要求解复杂非线性方程和方程组，求解过程是一个逼近过程。Matlab求解线性方程组 AX=B或XA=B 在MATLAB中，求解线性方程组时，主要采用前面章节介绍的除法运算符“/”和“\”。如： X=A\B表示求矩阵方程AX＝B的解； 的解。XA=B表示矩阵方程B/A＝X． 对方程组X＝A\B，要求A和B用相同的行数，X和B有相同的列数，它的行数等于矩阵A的列数，方程X＝B/A同理。 如果矩阵A不是方阵，其维数是m×n，则有： m＝n 恰定方程，求解精确解； m>n 超定方程，寻求最小二乘解； m 《MATLAB语言》课成论文 利用MATLAB求线性方程组 姓名：郭亚兰 学号：12010245331 专业：通信工程 班级：2010级通信工程一班 指导老师：汤全武 学院：物电学院 完成日期：2011年12月17日 利用MATLAB求解线性方程组 （郭亚兰 12010245331 2010 级通信一班） 【摘要】在高等数学及线性代数中涉及许多的数值问题，未知数的求解，微积分，不定积分，线性方程组的求解等对其手工求解都是比较复杂，而MATLAB语言正是处理线性方程组的求解的很好工具。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是讨论矩阵理论、与矩阵结合的有限维向量空间及其线性变换理论的一门学科。 【关键字】线性代数MATLAB语言秩矩阵解 一、基本概念 1、N级行列式A：A等于所有取自不同性不同列的n个元素的积的代数和。 2、矩阵B：矩阵的概念是很直观的，可以说是一张表。 3、线性无关：一向量组（a1,a2,…,an）不线性相关，既没有不全为零的数 k1,k2,………kn使得：k1*a1+k2*a2+………+kn*an=0 4、秩：向量组的极在线性无关组所含向量的个数成为这个向量组的秩。 5、矩阵B的秩：行秩，指矩阵的行向量组的秩；列秩类似。记：R(B) 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB 3-7变量非线性方程组的逆Broyden解法 matlab程序 function n_broyden clear all %清内存 clc %清屏 format long i=input('请输入非线性方程组个数(3-7)i='); switch i case 3 syms x y z x0=input('请输入初值（三维列向量[x;y;z])='); eps=input('请输入允许的误差精度eps='); f1=input('请输入f1(x,y,z)='); f2=input('请输入f2(x,y,z)='); f3=input('请输入f3(x,y,z)='); F=[f1;f2;f3]; J=jacobian(F,[x,y,z]); %求jacobi矩阵 Fx0=subs(F,{x,y,z},x0); Jx0=subs(J,{x,y,z},x0); H=inv(Jx0); x1=x0-H*Fx0; k=0; while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0; q=subs(F,{x,y,z},x1)-subs(F,{x,y,z},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1; x0=x1; Fx0=subs(F,{x,y,z},x0); x1=x1-H*Fx0; k=k+1; end x1 k case 4 syms a b c d x0=input('请输入初值（列向量[a;b;c;d])='); eps=input('请输入允许的误差精度eps='); f1=input('请输入f1(a,b,c,d)='); f2=input('请输入f2(a,b,c,d)='); f3=input('请输入f3(a,b,c,d)='); f4=input('请输入f4(a,b,c,d)='); F=[f1;f2;f3;f4]; J=jacobian(F,[a,b,c,d]); %求jacobi矩阵 Fx0=subs(F,{a,b,c,d},x0); Jx0=subs(J,{a,b,c,d},x0); H=inv(Jx0); x1=x0-H*Fx0; k=0; while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0; q=subs(F,{a,b,c,d},x1)-subs(F,{a,b,c,d},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1; x0=x1; Fx0=subs(F,{a,b,c,d},x0); x1=x1-H*Fx0; k=k+1; end x1 k case 5 syms a b c d e x0=input('请输入初值（列向量[a;b;c;d;e])='); eps=input('请输入允许的误差精度eps='); f1=input('请输入f1(a,b,c,d,e)='); f2=input('请输入f2(a,b,c,d,e)='); f3=input('请输入f3(a,b,c,d,e)='); f4=input('请输入f4(a,b,c,d,e)='); f5=input('请输入f5(a,b,c,d,e)='); F=[f1;f2;f3;f4;f5]; J=jacobian(F,[a,b,c,d,e]); %求jacobi矩阵 Fx0=subs(F,{a,b,c,d,e},x0); Jx0=subs(J,{a,b,c,d,e},x0); H=inv(Jx0); x1=x0-H*Fx0; 第7章 求解非线性方程 7.1 多项式运算在MATLAB 中的实现 一、多项式的表达 n 次多项式表达为：n a +??++=x a x a x a p(x )1-n 1-n 1n 0，是n+1项之和 在MATLAB 中，n 次多项式可以用n 次多项式系数构成的长度为n+1的行向量表示 [a0, a1,……an-1,an] 二、多项式的加减运算 设 有 两 个 多 项 式 n a +??++=x a x a x a p1(x )1-n 1-n 1n 0和 m b +??++=x b x b x b p2(x )1-m 1-m 1m 0。它们的加减运算实际上就是它们的对应系 数的加减运算。当它们的次数相同时，可以直接对多项式的系数向量进行加减运算。当它们的次数不同时，应该把次数低的多项式无高次项部分用0系数表示。 例2 计算()()1635223-+++-x x x x a=[1, -2, 5, 3]; b=[0, 0, 6, -1]; c=a+b 例3 设()6572532345++-+-=x x x x x x f ，()3532-+=x x x g ，求f(x)+g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; g1=[0, 0, 0, g];%为了和f 的次数找齐 f+g1, f-g1 三、多项式的乘法运算 conv(p1,p2) 例4 在上例中，求f(x)*g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; conv(f, g) 四、多项式的除法运算 [Q, r]=deconv(p1, p2) 表示p1除以p2，给出商式Q(x)，余式r(x)。Q,和r 仍为多项式系数向量 例4 在上例中，求f(x)/g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; [Q, r]=deconv(f, g) 五、多项式的导函数 p=polyder(P)：求多项式P 的导函数 p=polyder(P,Q)：求P·Q 的导函数 第4章 一阶线性微分方程组 一 内容提要 1． 基本概念 一阶微分方程组：形如 ??? ????? ???===) ,,,,( ),,,,(),,,,(2121222111 n n n n n y y y x f dx dy y y y x f dx dy y y y x f dx dy ΛΛΛΛΛ （3.1） 的方程组，（其中n y y y ,,,21Λ是关于x 的未知函数）叫做一阶微分方程组。 若存在一组函数)(,),(),(21x y x y x y n Λ使得在[a,b]上有恒等式 ),,2,1))((,),(),(,() (21n i x y x y x y x f dx x dy n i i ΛΛ==成立，则 )(,),(),(21x y x y x y n Λ称为一阶微分方程组(3.1)的一个解 含有n 任意常数n C C C ,,,21Λ的解 ?????? ?===) ,,,,( ),,,,(),,,,(21321222111n n n n C C C x y C C C x y C C C x y ΛΛΛΛΛ??? 称为(3.1)通解。如果通解满方程组 ???????=Φ=Φ=Φ0 ),,,,,,,,( 0),,,,,,,,(0),,,,,,,,(21212121221211n n n n n n n C C C y y y x C C C y y y x C C C y y y x ΛΛΛΛΛΛΛΛ 则称这个方程组为（3.1）的通积分。 满足初始条件,)(,,)(,)(0020021001n n y x y y x y y x y ===Λ的解，叫做初值问题的解。 实验1.1 用matlab 求解线性方程组 第一节 线性方程组的求解 一、齐次方程组的求解 rref （A ） %将矩阵A 化为阶梯形的最简式 null （A ） %求满足AX ＝0的解空间的一组基，即齐次线性方程组的基 础解系 【例１】 求下列齐次线性方程组的一个基础解系，并写出通解： 我们可以通过两种方法来解： 解法1： >> A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; >> rref(A) 执行后可得结果： ans= 1 -1 0 0 0 0 -1 1 0 0 0 0 由最简行阶梯型矩阵，得化简后的方程 ??? ??=+--=+--=-+-0 22004321 43214321x x x x x x x x x x x x 取x2，x4为自由未知量，扩充方程组为 即 提取自由未知量系数形成的列向量为基础解系，记 所以齐次方程组的通解为 解法2： clear A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; B=null(A, 'r') % help null 看看加个‘r’是什么作用, 若去掉r ，是什么结果？ 执行后可得结果： B= 1 0 1 0 0 1 0 1 ?? ?=-=-0 04321x x x x ?????? ?====4 4432221x x x x x x x x ??? ??? ??????+????????????=????? ???????1100001142 4321x x x x x x , 00111????? ? ??????=ε, 11002????? ???????=ε2 211εεk k x += 中南大学 MATLAB程序设计实践学长有爱奉献，下载填上信息即可上交，没有下载券 的自行百度。所需m文件照本文档做即可，即新建（FILE）→脚本（NEW-Sscript）→复制本文档代码→运行（会跳出 保存界面，文件名默认不要修改，保存）→结果。第 一题需要把数据文本文档和m文件放在一起。全部测 试无误，放心使用。本文档针对做牛顿法求非线性函 数题目的同学，当然第一题都一样，所有人都可以用。 ←记得删掉这段话 班级： 学号： 姓名： 一、《MATLAB程序设计实践》Matlab基础 表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散空间函数值来表示取向分布函数，是三维取向分布函数的一个实例。由于数据量非常大，不便于分析，需要借助图形来分析。请你编写一个matlab程序画出如下的几种图形来分析其取向分布特征：（1）用Slice函数给出其整体分布特征； （2）用pcolor或contour函数分别给出(φ2＝0, 5, 10, 15, 20, 25, 30, 35 … 90)切面上f分布情况(需要用到subplot函数)； (3) 用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。 备注：数据格式说明 解： （1）将文件内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如下： fid=fopen(''); for i=1:18 tline=fgetl(fid); end phi1=1;phi=1;phi2=1;line=0; f=zeros(19,19,19); while ~feof(fid) tline=fgetl(fid); data=str2num(tline); line=line+1; if mod(line,20)==1 phi2=(data/5)+1; phi=1; 数据说明部分，与作图无关 此方向表示f 随着φ1从0,5,10,15, 20 …到90的变化而变化 此方向表示f 随着φ从0,5,10,15, 20 …到90的变化而变化 表示以下数据为φ2＝0的数据，即f （φ1，φ，0） Matlab线性方程组求解 1. Gauss消元法： function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k); %计算行列式 end det=det*a(n,n); for k=n:-1:1 %回代求解 for j=k+1:n b(k)=b(k)-a(k,j)*x(j); end x(k)=b(k)/a(k,k); end Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 2. 列主元Gauss消去法： function x=detGauss(a,b) % Gauss列主元消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0; %选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n 基于Matlab 实现牛顿迭代法解非线性方程组 已知非线性方程组如下 2211221212 10801080x x x x x x x ?-++=??+-+=?? 给定初值0(0,0)T x =，要求求解精度达到0.00001 首先建立函数F(x)，方程组编程如下，将F.m 保存到工作路径中： function f=F(x) f(1)=x(1)^2-10*x(1)+x(2)^2+8; f(2)=x(1)*x(2)^2+x(1)-10*x(2)+8; f=[f(1) f(2)]; 建立函数DF(x),用于求方程组的Jacobi 矩阵，将DF.m 保存到工作路径中： function df=DF(x) df=[2*x(1)-10,2*x(2);x(2)^2+1,2*x(1)*x(2)-10]; 编程牛顿迭代法解非线性方程组，将newton.m 保存到工作路径中： clear; clc x=[0,0]'; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',0,x(1),x(2)); N=4; for i=1:N y=df\f'; x=x-y; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',i,x(1),x(2)); if norm(y)<0.0000001 break ; else end end 运行结果如下： 0 0.0000000 0.0000000 1 0.8000000 0.8800000 2 0.9917872 0.9917117 3 0.9999752 0.9999685 4 1.0000000 1.0000000 遗传算法解非线性方程组的Matlab程序 程序用MATLAB语言编写。之所以选择MATLB，是因为它简单，但又功能强大。写1行MATLAB程序，相当于写10行C++程序。在编写算法阶段，最好用MATLAB语言，算法验证以后，要进入工程阶段，再把它翻译成C++语言。 本程序的算法很简单，只具有示意性，不能用于实战。 非线性方程组的实例在函数(2)nonLinearSumError1(x)中，你可以用这个实例做样子构造你自己待解的非线性方程组。 %注意：标准遗传算法的一个重要概念是，染色体是可能解的2进制顺序号，由这个序号在可能解的集合(解空间)中找到可能解 %程序的流程如下： %程序初始化，随机生成一组可能解(第一批染色体) %1: 由可能解的序号寻找解本身(关键步骤) %2：把解代入非线性方程计算误差，如果误差符合要求，停止计算 %3：选择最好解对应的最优染色体 %4：保留每次迭代产生的最好的染色体，以防最好染色体丢失 %5: 把保留的最好的染色体holdBestChromosome加入到染色体群中 %6: 为每一条染色体(即可能解的序号)定义一个概率(关键步骤) %7：按照概率筛选染色体(关键步骤) %8：染色体杂交(关键步骤) %9：变异 %10:到1 %这是遗传算法的主程序，它需要调用的函数如下。 %由染色体(可能解的2进制)顺序号找到可能解： %(1)x=chromosome_x(fatherChromosomeGroup,oneDimensionSet,solutionSum); %把解代入非线性方程组计算误差函数：(2)functionError=nonLinearSumError1(x); %判定程是否得解函数：(3)[solution,isTrue]=isSolution(x,funtionError,solutionSumError); %选择最优染色体函数： %(4)[bestChromosome,leastFunctionError]=best_worstChromosome(fatherChromosomeGroup,functionError); %误差比较函数：从两个染色体中，选出误差较小的染色体 %(5)[holdBestChromosome,holdLeastFunctionError]... % =compareBestChromosome(holdBestChromosome,holdLeastFunctionError,... % bestChromosome,leastFuntionError) %为染色体定义概率函数，好的染色体概率高，坏染色体概率低 %(6)p=chromosomeProbability(functionError); %按概率选择染色体函数： %(7)slecteChromosomeGroup=selecteChromome(fatherChromosomeGroup,p); %父代染色体杂交产生子代染色体函数 %(8)sonChrmosomeGroup=crossChromosome(slecteChromosomeGroup,2); %防止染色体超出解空间的函数 %(9)chromosomeGroup=checkSequence(chromosomeGroup,solutionSum) %变异函数 %(10)fatherChromosomeGroup=varianceCh(sonChromosomeGroup,0.8,solutionN); %通过实验有如下结果： %1。染色体应当多一些 本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（M üller ）法和迭代法的加速等）及其MATLAB 程序，求解非线性方程组的方法及其MATLAB 程序. 2.1 方程（组）的根及其MATLAB 命令 2.1.2 求解方程（组）的solve 命令 求方程f (x )=q (x )的根可以用MATLAB 命令: >> x=solve('方程f(x)=q(x)',’待求符号变量x ’) 求方程组f i (x 1,…,x n )=q i (x 1,…,x n ) (i =1,2,…,n )的根可以用MATLAB 命令: >>E1=sym('方程f1(x1,…,xn)=q1(x1,…,xn)'); ……………………………………………………. En=sym('方程fn(x1,…,xn)=qn(x1,…,xn)'); [x1,x2,…,xn]=solve(E1,E2,…,En, x1,…,xn) 2.1.3 求解多项式方程（组）的roots 命令 如果)(x f 为多项式，则可分别用如下命令求方程0)(=x f 的根，或求导数)('x f （见表 2-1）. 2.1.4 求解方程（组）的fsolve 命令 如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots 命令.如果非线性方程（组）是含有超越函数，则无法使用roots 命令，需要调用MATLAB 系统中提供的另一个程序fsolve 来求解.当然，程序fsolve 也可以用于多项式方程（组），但是它的计算量明显比roots 命令的大. fsolve 命令使用最小二乘法（least squares method ）解非线性方程（组） (F X =)0 的数值解，其中X 和F （X ）可以是向量或矩阵.此种方法需要尝试着输入解X 的初始值（向量或矩阵）X 0，即使程序中的迭代序列收敛，也不一定收敛到(F X =)0的根（见例2.1.8）. fsolve 的调用格式： X=fsolve(F,X0) 输入函数)(x F 的M 文件名和解X 的初始值（向量或矩阵）X 0，尝试着解方程（组） 第四讲 常系数线性微分方程组的解法（4课时） 一、目的与要求: 理解常系数线性微分方程组的特征方程式, 特征根, 特征向量的概念, 掌 握常系数线性微分方程组的基本解组的求法. 二、重点：常系数线性微分方程组的基本解组的求法. 三、难点：常系数线性微分方程组的特征方程式, 特征根, 特征向量的概念. 四、教学方法：讲练结合法、启发式与提问式相结合教学法. 五、教学手段：传统板书与多媒体课件辅助教学相结合. 六、教学过程： 1 新课引入 由定理3.6我们已知道，求线性齐次方程组(3.8)的通解问题，归结到求其基本解组. 但是对于一般的方程组(3.8)，如何求出基本解组，至今尚无一般方法. 然而对于常系数线性齐次方程组 dY AY dx = (3.20) 其中A 是n n ?实常数矩阵，借助于线性代数中的约当(Jordan)标准型理论或矩阵指数，可以使这一问题得到彻底解决. 本节将介绍前一种方法，因为它比较直观. 由线性代数知识可知，对于任一n n ?矩阵A ，恒存在非奇异的n n ?矩阵T ，使矩阵 1T AT -成为约当标准型. 为此，对方程组(3.20)引入非奇异线性变换 Y TZ = (3.21) 其中()(,1,2,,),ij T t i j n == det 0T ≠,将方程组(3.20)化为 1dZ T ATZ dx -= (3.22) 我们知道，约当标准型1 T AT -的形式与矩阵A 的特征方程 11121212221 2 det()0n n n n nn a a a a a a A E a a a λ λλλ ---= =- 的根的情况有关. 上述方程也称为常系数齐次方程组(3.20)的特征方程式.它的根称为矩阵 A 的特征根. 下面分两种情况讨论. (一) 矩阵A 的特征根均是单根的情形. 设特征根为12,,,,n λλλ 这时 12 1 00 n T AT λλλ-????? ?=?????? 方程组(3.20)变为 11122 200n n n dz dx z dz z dx z dz dx λλλ?????????????? ????????= ???????????????? ?????? (3.23) 易见方程组(3.23)有n 个解 1110(),00x Z x e λ????????=???????? 220010(),,()0001n x x n Z x e Z x e λλ???????????? ????==???????????????? 把这n 个解代回变换(3.21)之中，便得到方程组(3.20)的n 个解 12()i i i i x x i i ni t t Y x e e T t λλ?? ????==?????? (1,2,,)i n = 实用数值方法(Matlab) 综述报告题目：非线性方程（组）求根问题 小组成员 许多数学和物理问题归结为解函数方程f(x)=0。方程f(x)=0的解称为方程的根。对于非线性方程，在某个范围内往往不止一个根，而且根的分布情况可能很复杂，面对这种情况，通常先将考察的范围花费为若干个子段，然后判断哪些子段内有根，然后再在有根子段内找出满足精度要求的近似根。为此适当选取有根子段内某一点作为根的初始值近似，然后运用迭代方法使之足部精确化。这就是方程求根的迭代法。下面介绍书上的几种方法： 1、二分法 （1）方法概要： 假定函数f(x)在[a,b]上连续，且f(a)f(b)=0,则方程f(x)=0在[a,b]内一定有实根。取其中 将其二分，判断所求的根在的左侧还是右侧，得到一个新的有根区间 点 [],长度为[a,b]的一半。对新的有根区间继续实行上述二分手段，直至二分k次后有根区间[]长度 可见，如果二分过程无限继续下去，这些有限根区间最终必收敛于一点,该点就是所求的根。在实际计算过程中不可能完成这个无限过程，允许有一定的误差，则二分k+1次后 只要有根区间[]的长度小于，那么结果关于允许误差就能“准确”地满足方程f(x)=0。 (2)计算框图： 2、开方法 对于给定,求开方值 为此，可以运用校正技术设计从预报值生成校正值的迭代公式。自然希望校正值 能更好满足所给方程： 这是个关于校正量的近似关系式，如果从中删去二次项，即可化归为一次方程 解之有 从而关于校正值有如下开方公式 上述演绎过程表明，开方法的设计思想是逐步线性化，即将二次方程的求解画归为一次方程求解过程的重复。开方公式规定了预报值与校正值之间的一种函数关系,这里 为开方法的迭代函数。 3、Newton法 （1）方法概要 线性方程组求解 1.直接法 Gauss消元法： function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k); end det=det*a(n,n); for k=n:-1:1 %回代 for j=k+1:n b(k)=b(k)-a(k,j)*x(j); end x(k)=b(k)/a(k,k); end Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 列主元Gauss消去法： function x=detGauss(a,b) % Gauss列主元消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0;% 选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n z=a(k,j);a(k,j)=a(r,j);a(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det; endMATLAB解线性方程组的直接方法
Matlab求解线性方程组非线性方程组
利用MATLAB求线性方程组
线性方程组求解matlab实现
3-7变量非线性方程组的逆Broyden解法matlab程序
MATLAB应用 求解非线性方程
一阶线性微分方程组
实验一用matlab求解线性方程组
matlab程序设计实践-牛顿法解非线性方程
Matlab线性方程组求解(Gauss消去法)
基于Matlab的牛顿迭代法解非线性方程组
遗传算法解非线性方程组的Matlab程序
第二章非线性方程(组)的数值解法的matlab程序
第三章 一阶线性微分方程组 第四讲 常系数线性微分方程组的解法(1)
MATLAB 非线性方程(组)求根
线性方程组求解Matlab程序(精.选)
相关主题
文本预览