下载文档原格式
实验二 卷积及微分、差分方程的求解一、实验目的1．学习利用Matlab 工程软件实现系统时域输入—输出分析法2．观察和掌握各种常用信号的波形及其卷积后的结果。
3．利用Matlab 实现微分、差分方程的求解。
二、原理说明在连续线性时不变因果（LTI ）系统的时域分析中，我们可以通过经典求解微分方程的方法，同样在离散线性时不变因果（LTI ）系统的时域分析中我们也可以利用经典的差分方程求解法。
这两种方法都是通过求解方程的齐次解和特解得到的，无论哪种方法都需要经过精确、繁琐的数学计算，它们在高等数学课中都有详细介绍，这里我们不再赘述。
对于线性系统的零状态响应r f (t)，我们定义为在初始状态为零、系统只在激励信号e(t)作用下所产生的输出，若该系统的单位冲激响应为h(t)[系统在单位冲激信号作用下的零状态响应]，则系统的⎰∞--∙=*=tf d t h e t h t e t r τττ)()()()()( 引入Matlab 算法以后，我们可以利用有限的几个函数就可以将方程的输入—输出关系求出，并且以图形的方式表达，从而方便、直接地观察系统的时域特性。
三、预习要求1．所用函数（1）单位阶跃信号：Heaviside(t)；功能：产生单位阶跃信号。
（2）门函数：rectpuls(t);功能：产生门函数。
（3）函数square功能：产生周期为2幅值为1的方波信号，调用格式：x=square(t) x=square(t,duty)，其中，t 为时间向量duty 为正幅值部分占周期的百分比。
（4）sawtooth ：功能：产生锯齿波或三角波；调用格式：x=sawtooth(t) x=sawtooth(t,width)；sawtooth(t)用于产生周期为2π，最大幅度为1，最小幅度为-1的周期性三角波（锯齿波），其中，width 参数为0和1之间的数，表示最大幅度出现的位置。
Width=0.5时，产生标准正三角波。
实验二 卷积及微分、差分方程的求解一、实验目的1．学习利用Matlab 工程软件实现系统时域输入—输出分析法2．观察和掌握各种常用信号的波形及其卷积后的结果。
3．利用Matlab 实现微分、差分方程的求解。
二、原理说明在连续线性时不变因果（LTI ）系统的时域分析中，我们可以通过经典求解微分方程的方法，同样在离散线性时不变因果（LTI ）系统的时域分析中我们也可以利用经典的差分方程求解法。
这两种方法都是通过求解方程的齐次解和特解得到的，无论哪种方法都需要经过精确、繁琐的数学计算，它们在高等数学课中都有详细介绍，这里我们不再赘述。
对于线性系统的零状态响应r f (t)，我们定义为在初始状态为零、系统只在激励信号e(t)作用下所产生的输出，若该系统的单位冲激响应为h(t)[系统在单位冲激信号作用下的零状态响应]，则系统的⎰∞--∙=*=tf d t h e t h t e t r τττ)()()()()( 引入Matlab 算法以后，我们可以利用有限的几个函数就可以将方程的输入—输出关系求出，并且以图形的方式表达，从而方便、直接地观察系统的时域特性。
三、预习要求1．所用函数（1）单位阶跃信号：Heaviside(t)；功能：产生单位阶跃信号。
（2）门函数：rectpuls(t);功能：产生门函数。
（3）函数square功能：产生周期为2幅值为1的方波信号，调用格式：x=square(t) x=square(t,duty)，其中，t 为时间向量duty 为正幅值部分占周期的百分比。
（4）sawtooth ：功能：产生锯齿波或三角波；调用格式：x=sawtooth(t) x=sawtooth(t,width)；sawtooth(t)用于产生周期为2π，最大幅度为1，最小幅度为-1的周期性三角波（锯齿波），其中，width 参数为0和1之间的数，表示最大幅度出现的位置。
Width=0.5时，产生标准正三角波。
微分方程的数值解法微分方程的数值解法微分方程的数值解法【1】摘要：本文结合数例详细阐述了最基本的解决常微分方程初值问题的数值法，即Euler方法、改进Euler法，并进行了对比，总结了它们各自的优点和缺点，为我们深入探究微分方程的其他解法打下了坚实的基础。
关键词:常微分方程数值解法 Euler方法改进Euler法1、Euler方法由微分方程的相关概念可知，初值问题的解就是一条过点的积分曲线，并且在该曲线上任一点处的切线斜率等于函数的值。
根据数值解法的基本思想，我们取等距节点，其中h为步长，在点处，以为斜率作直线交直线于点。
如果步长比较小，那么所作直线与曲线的偏差不会太大，所以可用的近似值，即：，再从点出发，以为斜率作直线，作为的近似值，即：重复上述步骤，就能逐步求出准确解在各节点处的近似值。
一般地，若为的近似值，则过点以为斜率的直线为：从而的近似值为：此公式就是Euler公式。
因为Euler方法的思想是用折线近似代替曲线，所以Euler方法又称Euler折线法。
Euler方法是初值问题数值解中最简单的一种方法，由于它的精度不高，当步数增多时，由于误差的积累，用Euler方法作出的折线可能会越来越偏离曲线。
举例说明：解： ,精确解为：1.2 -0.96 -1 0.041.4 -0.84 -0.933 0.9331.6 -0.64 -0.8 0.161.8 -0.36 -0.6 0.242.0 0 -0.333 0.332.2 0.44 0 0.44通过上表可以比较明显地看出误差随着计算在积累。
2、改进Euler法方法构造在常微分方程初值问题 ,对其从到进行定积分得:用梯形公式将右端的定积分进行近似计算得：用和来分别代替和得计算格式:这就是改进的Euler法。
解：解得：由于 ,是线形函数可以从隐式格式中解出问题的精确解是误差0.2 2.421403 2.422222 0.000813 0.021400.4 2.891825 2.893827 0.00200 0.051830.6 3.422119 3.425789 0.00367 0.094112.0 10.38906 10.43878 0.04872 1.1973通过比较上表的第四列与第五列就能非常明显看出改进Euler方法精度比Euler方法精度高。
讲义二：卷积与微分方程的数值法求解一、 从离散卷积和到连续卷积序列f 1(k )和f 2(k )的离散卷积定义式为()()()()1212i f k f k f i f k i ∞=−∞∗=−∑ 用来计算离散卷积的函数为：f=conv(f1,f2) f1,f2为参与卷积运算的两个序列，f 为卷积的结果，长度为length(f1)+length(f2)-1。
[f,r]=deconv(f1,f2) 解卷运算，使f1=conv(f,f2)+r 成立EX 错误！文档中没有指定样式的文字。
-1 ()()1sin ,010x k k k =≤≤，()20.8,015k x k k =≤≤，计算离散卷积和()y k =()1x k ∗()2x k 。
%程序5_1 计算离散卷积和k1=0:10; %x1的变量取值范围x1=sin(k1); %构建x1序列k2=0:15; %x2的变量取值范围x2=0.8.^k2; %构建x2序列y=conv(x1,x2); %计算卷积结果%显示卷积结果subplot(3,1,1);stem(k1,x1);title('x_1(k)');subplot(3,1,2);stem(k2,x2);title('x_2(k)');k=0:length(y)-1;subplot(3,1,3);stem(k,y);title('y(k)');下面讨论连续卷积的计算：连续时间函数1()f t 和2()f t 的卷积定义为：()()()()()1212f t f t f t f f t d τττ∞−∞=∗=−∫由于计算机实际处理的数据必须满足：1、离散存储；2、有限数据量。
连续信号的处理必须首先经过数值化的过程，以离散的形式被分析、保存和处理。
用数值方法计算卷积需要将卷积积分看作信号的分段求和来实现，这样会得到一定的精确度要求下的卷积。
()()()()()()()1212120lim k f t f t f t f f t d f k f t k τττ∞∞−∞Δ→=−∞=∗=−=Δ−ΔΔ∑∫ 如果我们只求当t n =Δ(n 为整数)时f (t )的值()f n Δ,则得：()()()()1212[()]k k f n f k f n k f k f n k ∞∞=−∞=−∞Δ≈ΔΔ−ΔΔ=ΔΔ−Δ∑∑ 式中的()12[()]k f k f n k ∞=−∞Δ−Δ∑实际上就是连续信号f 1(t )和f 2(t )经等时间间隔Δ均匀抽样的离散序列1()f k Δ和2()f k Δ的离散卷积和。
微分方程的数值解法微分方程（Differential Equation）是描述自然界中变化的现象的重要工具，具有广泛的应用范围。
对于一般的微分方程，往往很难找到解析解，这时候就需要使用数值解法来近似求解微分方程。
本文将介绍几种常见的微分方程数值解法及其原理。
一、欧拉方法（Euler's Method）欧拉方法是最基本也是最容易理解的数值解法之一。
它的基本思想是将微分方程转化为差分方程，通过给定的初始条件，在离散的点上逐步计算出函数的近似值。
对于一阶常微分方程dy/dx = f(x, y)，利用欧拉方法可以得到近似解：y_n+1 = y_n + h * f(x_n, y_n)其中，h是步长，x_n和y_n是已知点的坐标。
欧拉方法的优点在于简单易懂，但是由于是一阶方法，误差较大，对于复杂的微分方程可能不够准确。
二、改进的欧拉方法（Improved Euler's Method）改进的欧拉方法又称为改进的欧拉-柯西方法，是对欧拉方法的一种改进。
它通过在每一步计算中利用两个不同点的斜率来更准确地逼近函数的值。
对于一阶常微分方程dy/dx = f(x, y)，改进的欧拉方法的迭代公式为：y_n+1 = y_n + (h/2) * [f(x_n, y_n) + f(x_n+1, y_n + h * f(x_n, y_n))]相较于欧拉方法，改进的欧拉方法具有更高的精度，在同样的步长下得到的结果更接近真实解。
三、四阶龙格-库塔方法（Fourth-Order Runge-Kutta Method）四阶龙格-库塔方法是一种更高阶的数值解法，通过计算多个点的斜率进行加权平均，得到更为准确的解。
对于一阶常微分方程dy/dx = f(x, y)，四阶龙格-库塔方法的迭代公式为：k1 = h * f(x_n, y_n)k2 = h * f(x_n + h/2, y_n + k1/2)k3 = h * f(x_n + h/2, y_n + k2/2)k4 = h * f(x_n + h, y_n + k3)y_n+1 = y_n + (k1 + 2k2 + 2k3 + k4)/6四阶龙格-库塔方法是数值解法中精度最高的方法之一，它的计算复杂度较高，但是能够提供更为准确的结果。
微分方程是数学中的一种重要的方程类型，广泛应用于物理、工程、经济等领域。
解微分方程有各种方法，其中数值解法是一种重要而实用的方法。
微分方程的数值解法是通过数值计算来求解微分方程的近似解。
它的基本思想是将微分方程转化为差分方程，并用计算机进行迭代计算，从而求得微分方程的数值解。
数值解法的关键在于如何将微分方程转化为差分方程。
常见的方法有欧拉方法、改进欧拉方法、龙格－库塔方法等。
这些方法都是基于泰勒级数展开的原理进行推导的。
以欧拉方法为例，其基本思路是将微分方程中的导数用差商的方式近似表示，然后通过迭代计算，逐步逼近微分方程的解。
欧拉方法的具体步骤如下：首先确定微分方程的初始条件，即给定t0时刻的函数值y0，然后选取一定的步长ℎ，利用微分方程的导数计算差商y′=dy，进而根据差商dt得到下一个时刻的函数值y n+1=y n+ℎy′。
通过不断迭代计算，即可得到微分方程在一定时间区间内的数值解。
数值解法的另一个重要问题是误差控制。
由于数值计算本身的误差以及近似方法的误差，数值解法所得到的结果通常与真实解存在误差。
为了控制误差，常用的方法有缩小步长ℎ、提高近似方法的阶数等。
此外，还可以通过与解析解进行比较，评估数值解的准确性。
微分方程的数值解法具有以下几点优势。
首先，微分方程的解析解通常较难求得，而数值解法可以给出一个近似解，提供了一种有效的解决方案。
其次，数值解法可以利用计算机的高速运算能力，进行大规模复杂微分方程的求解。
此外，数值解法还可以在实际问题中进行仿真和优化，即通过调整参数来求解微分方程，从而得到最优解。
尽管微分方程的数值解法具有广泛的应用前景，但也存在一些问题和挑战。
首先，数值解法的稳定性和收敛性需要深入研究和分析。
其次，数值解法的计算量通常较大，对计算机运算能力和存储空间的要求较高。
此外，数值解法还需要对问题进行适当的离散化处理，从而可能引入一定的误差。
综上所述，“微分方程的数值解法”是一种重要而实用的方法，可以有效地求解微分方程的近似解。
微分方程数值解法微分方程是数学中的重要概念，它描述了物理系统中变量之间的关系。
解微分方程是许多科学领域中常见的问题，其中又可以分为解析解和数值解两种方法。
本文将重点介绍微分方程的数值解法，并详细讨论其中的常用方法和应用。
一、微分方程的数值解法概述微分方程的解析解往往较为复杂，难以直接求解。
在实际问题中，我们通常利用计算机进行数值计算，以获得方程的数值解。
数值解法的基本思想是将微分方程转化为一组离散的数值问题，通过逼近连续函数来获得数值解。
二、常见的数值解法1. 欧拉法欧拉法是最基础的数值解法之一，其核心思想是将微分方程转化为差分方程，通过逼近连续函数来获得数值解。
欧拉法的基本形式为：yn+1 = yn + h·f(xn, yn)其中，yn表示第n个时间步的数值解，h为时间步长，f为微分方程右端的函数。
欧拉法的精度较低，但计算简单，适用于初步估计或简单系统的求解。
2. 改进的欧拉法（Heun法）改进的欧拉法（Heun法）是对欧拉法的改进，其关键在于求解下一个时间步的近似值时，利用了两个斜率的平均值。
Heun法的基本形式为：yn+1 = yn + (h/2)·(k1 + k2)k1 = f(xn, yn)，k2 = f(xn+h, yn+h·k1)Heun法较欧拉法的精度更高，但计算量较大。
3. 龙格-库塔法（RK方法）龙格-库塔法是一类常用的数值解法，包含了多个不同阶数的方法。
其中，最常用的是经典四阶龙格-库塔法（RK4法），其基本形式为：k1 = f(xn, yn)k2 = f(xn + h/2, yn + (h/2)·k1)k3 = f(xn + h/2, yn + (h/2)·k2)k4 = f(xn + h, yn + h·k3)yn+1 = yn + (h/6)·(k1 + 2k2 + 2k3 + k4)RK4法实现较为复杂，但精度较高，适用于解决大多数常微分方程问题。
讲义二：卷积与微分方程的数值法求解
一、 从离散卷积和到连续卷积
序列f 1(k )和f 2(k )的离散卷积定义式为
()()()()1212i f k f k f i f k i ∞=−∞∗=
−∑ 用来计算离散卷积的函数为：
f=conv(f1,f2) f1,f2为参与卷积运算的两个序列，f 为卷积的
结果，长度为length(f1)+length(f2)-1。
[f,r]=deconv(f1,f2) 解卷运算，使f1=conv(f,f2)+r 成立
EX 错误！文档中没有指定样式的文字。
-1 ()()1sin ,010x k k k =≤≤，
()20.8,015k x k k =≤≤，计算离散卷积和()y k =()1x k ∗()2x k 。
%程序5_1 计算离散卷积和
k1=0:10; %x1的变量取值范围
x1=sin(k1); %构建x1序列
k2=0:15; %x2的变量取值范围
x2=0.8.^k2; %构建x2序列
y=conv(x1,x2); %计算卷积结果
%显示卷积结果
subplot(3,1,1);stem(k1,x1);title('x_1(k)');
subplot(3,1,2);stem(k2,x2);title('x_2(k)');
k=0:length(y)-1;
subplot(3,1,3);stem(k,y);title('y(k)');
下面讨论连续卷积的计算：
连续时间函数1()f t 和2()f t 的卷积定义为：
()()()()()1212f t f t f t f f t d τττ∞
−∞=∗=−∫
由于计算机实际处理的数据必须满足：1、离散存储；2、有限数据量。
连续信号的处理必须首先经过数值化的过程，以离散的形式被分析、保存和处理。
用数值方法计算卷积需要将卷积积分看作信号的分段求和来实现，这样会得到一定的精确度要求下的卷积。
()()()()()()()1212120lim k f t f t f t f f t d f k f t k τττ∞
∞
−∞Δ→=−∞=∗=−=Δ−ΔΔ∑∫ 如果我们只求当t n =Δ(n 为整数)时f (t )的值()f n Δ,则得：
()()()()1212[()]k k f n f k f n k f k f n k ∞∞=−∞=−∞Δ≈
ΔΔ−ΔΔ=ΔΔ−Δ∑∑ 式中的()12[()]k f k f n k ∞=−∞
Δ−Δ∑实际上就是连续信号f 1(t )和f 2(t )经等时间
间隔Δ均匀抽样的离散序列1()f k Δ和2()f k Δ的离散卷积和。
当Δ足够小
时，()f n Δ就是卷积积分的近似计算结果。
因此，MATLAB 实现数值卷积的过程可总结如下：
1． 将连续信号f 1(t )和f 2(t )以时间间隔Δ进行取样，得到离散序列1()f k Δ和2()f k Δ。
2． 构造与1()f k Δ和2()f k Δ相对应的时间向量k1 和k2(向量k1和k2 是
取样时间间隔 的整数倍)。
3． 调用MATLAB 计算离散卷积和的指令conv()函数计算卷积积分f (t )的近似向量()f n Δ。
4． 构造()f n Δ对应的时间向量k(共length(k1)+length(k2)-1项，由k1与k2首项之和开始，以Δ为间隔，至k1与k2长度之和减2为止)。
例：数值计算f 1(t )=e -t ε(t )，f 2(t )=te -t ε(t )的卷积。
%程序2-6 卷积的数值计算
dt=0.01;
k1=0:dt:6;
f1=exp(-k1);%生成信号f1
k2=k1;
f2=k2.*exp(-k2);%生成信号f2
f=dt*conv(f1,f2);%计算卷积结果f
k0=k1(1)+k2(1); %计算序列f 非零样值的起点位置
k3=length(f1)+length(f2)-2; %计算卷积和f 的非零样值的宽度
k=k0:dt:k0+k3*dt; %确定卷积和f 非零样值的时间向量 subplot(2,2,1);
plot(k1,f1);title('f1(t)');xlabel('t'); %在子图1 绘f1(t)时域波形图 subplot(2,2,2);
plot(k2,f2);title('f2(t)');xlabel('t'); %在子图2 绘f2(t)时波形图 subplot(2,2,3);
plot(k,f); %画卷积f(t)的时域波形
h=get(gca,'position');
h(3)=2.5*h(3);set(gca,'position',h); %将第三个子图的横坐标范围扩为原来的2.5 倍
title('f(t)=f1(t)*f2(t)');xlabel('t');
二、 微分方程的数值求解
基本思想是把微分方程化为差分方程即可。
比如二阶微分方程：''''y ay by cf df ++=+，两个初始条件一般设定为：y(0-)；y'(0-)。
利用计算机进行数值化近似求解时，需要把微分化为差分：
()()y t y k ⇔
()1[()(1)]s
dy t y k y k dt T ⇔−− 222
()1[()2(1)(2)]s d y t y k y k y k dt T ⇔−−+− 利用这些转化公式，原微分方程就变成了一个形如：
()(1)(2)()(1)y k y k y k f k f k αβγλ+−+−=+−形式的后向差分方程。
求解原微分方程就转化为求解差分方程。
但是求解差分方程需要的初始条件y(-
1)；y(-2)如何得到呢？
考虑到y(0-)是激励加入前连续系统的初始状态，相应离散系统激励加入前的初始状态是y(-1)，所以又可以得到以下方程组： (0)(1)
(1)(2)'(0)S
y y y y y T −≈−−−−−≈ 根据此式求出差分方程的初始条件，之后用迭代法求解即可。
如果连续系统给出的是初始条件y(0+)；y'(0+)，那么求离散系统初始条件的方程组应该为：
(0)(0)
(1)(0)'(0)S
y y y y y T +≈−+≈
例：数值法求解微分方程''4'3'3y y y f f ++=+，初始条件为：y(0-)=y'(0-)=1，()()f t t ε=。
步长取T S =0.1，求解范围为区间[0,4]。
解：计算知全响应为：3()(1)()t t y t e e t ε−−=−+，绘制曲线进行对照。
把微分方程转化为差分方程为
143()240(1)100(2)13()10(1)y k y k y k f k f k −−+−=−−
初始条件需要求解方程组：
(1)1
(1)(2)10.1
y y y −=−−−= 得到：y(-1)=1；y(-2)=0.9。
Matlab 程序如下．
Ts=0.1; %时间间隔
n=1:4/Ts+2; %坐标1,2代表初始条件坐标-2,-1
F(n)=1;F(1)=0;F(2)=0; %激励信号初始条件
Y(1)=0.9;Y(2)=1; %输入初始条件y(-2)=0.9;y(-1)=1
for k=3:length(n) %迭代求解
Y(k)=(13*F(k)-10*F(k-1)-100*Y(k-2)+240*Y(k-1))/143;
end
n1=(n-3)*Ts;
subplot(2,1,1);plot(n1,Y);
xlabel('t');ylabel('y(t)');title('全响应(数值计算)');grid on;
subplot(2,1,2);ezplot('-exp(-3*t)+exp(-t)+1',[0,4]);
xlabel('t');ylabel('y(t)');title('全响应');grid on;
为了使用方便，MATLAB 提供了直接数值解算微分方程的常用指令如：ode45、ode23等；利用符号工具包中的dsolve 函数也可以很方便地求解微分方程。
这些指令求解的微分方程不限于线性常系数微分方程。
ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：
ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．
ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度． 例：求解微分方程初值问题2222dy y x x dx
+=+，(0)1y −=的数值解，求解范围为区间[0, 0.5]．
fun=inline('-2*y+2*x^2+2*x','x','y');
[x,y]=ode23(fun,[0,0.5],1);
x';
y';
plot(x,y,'o-')。
