差分方程及matlab求解

合集下载

怎样用Matlab求解差分方程题解读

模型建立记第k年沙丘鹤的数量为xk,年均增长率为 r，则第k+1年鹤的数量为
xk+1=(1+r)xk k=0,1,2······
已知x0=100, 在较好，中等和较差的自然环境下 r=0.0194, -0.0324,和-0.0382 我们利用 Matlab编程，递推20年后观察沙丘鹤的数量变化情况
Xk-1决定的部分是 a1bcXk-1,由Xk-2决定的部分是 a2b(1-a1)bcXk-2
Xk= a1bcXk-1 + a2b(1-a1)bcXk-2
Xk= a1bcXk-1 + a2b(1-a1)bcXk-2
实际上，就是Xk= pXk-1 + qXk-2 我们需要知道x0,a1,a2,c, 考察b不同时，种子繁殖的情况。在这里假设 X0=100,a1=0.5,a2=0.25,c=10,b=0.18~0.20 这样可以用matlab计算了
K=(0:20)’; Y1=zwfz(100,21,0.18); Y2=zwfz(100,21,0.19); Y3=zwfz(100,21,0,20); Round([k,y1’,y2’,y3’]) Plot(k,y1,k,y2,’:’,k,y3,’o’), Gtext(‘b=0.18’),gtext(‘b=0.19’),gtext(‘b=0.20’)
function x=czqc(n) A=[0.6,0.2,0.1;0.3,0.7,0.3;0.1,0.1,0.6]; x(:,1)=[200,200,200]'; for k=1:n x(:,k+1)=A*x(:,k); end
如果直接看10年或者20年发展趋势，可以直接在命令窗口（commond window）作，而不是必须编一个函数

实验二微分方程与差分方程模型Matlab求解

实验二：微分方程与差分方程模型Matlab 求解一、实验目的[1] 掌握解析、数值解法，并学会用图形观察解的形态和进展解的定性分析；[2] 熟悉MATLAB 软件关于微分方程求解的各种命令； [3] 通过范例学习建立微分方程方面的数学模型以及求解全过程；[4] 熟悉离散 Logistic 模型的求解与混沌的产生过程。

二、实验原理1. 微分方程模型与MATLAB 求解解析解用MATLAB 命令dsolve(‘eqn1’,’eqn2’, ...) 求常微分方程〔组〕的解析解。

其中‘eqni'表示第i 个微分方程，Dny 表示y 的n 阶导数,默认的自变量为t 。

〔1〕微分方程例1 求解一阶微分方程 21y dxdy+= (1) 求通解输入：dsolve('Dy=1+y^2') 输出： ans =tan(t+C1) 〔2〕求特解输入：dsolve('Dy=1+y^2','y(0)=1','x') 指定初值为1，自变量为x输出： ans =tan(x+1/4*pi)例2 求解二阶微分方程 221()04(/2)2(/2)2/x y xy x y y y πππ'''++-=='=-原方程两边都除以2x ，得211(1)04y y y x x '''++-= 输入:dsolve('D2y+(1/x)*Dy+(1-1/4/x^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x') ans =- (exp(x*i)*(pi/2)^(1/2)*i)/x^(1/2) +(exp(x*i)*exp(-x*2*i)*(pi/2)^(3/2)*2*i)/(pi*x^(1/2))试试能不用用simplify 函数化简输入: simplify(ans)ans =2^(1/2)*pi^(1/2)/x^(1/2)*sin(x) 〔2〕微分方程组例3 求解 d f /d x =3f +4g ; d g /d x =-4f +3g 。

差分方程的解法分析及MATLAB实现(程序)

差分方程的解法分析及MATLAB 实现（程序）摘自：张登奇,彭仕玉.差分方程的解法分析及其MATLAB 实现[J]. 湖南理工学院学报.2014(03) 引言线性常系数差分方程是描述线性时不变离散时间系统的数学模型,求解差分方程是分析离散时间系统的重要内容.在《信号与系统》课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域法[1].1 迭代法例1 已知离散系统的差分方程为)1(31)()2(81)1(43)(-+=-+--n x n x n y n y n y ,激励信号为)()43()(n u n x n =,初始状态为21)2(4)1(=-=-y y ，.求系统响应. 根据激励信号和初始状态,手工依次迭代可算出2459)1(,25)0(==y y . 利用MATLAB 中的filter 函数实现迭代过程的m 程序如下:clc;clear;format compact;a=[1,-3/4,1/8],b=[1,1/3,0], %输入差分方程系数向量,不足补0对齐n=0:10;xn=(3/4).^n, %输入激励信号zx=[0,0],zy=[4,12], %输入初始状态zi=filtic(b,a,zy,zx),%计算等效初始条件[yn,zf]=filter(b,a,xn,zi),%迭代计算输出和后段等效初始条件2 时域经典法用时域经典法求解差分方程:先求齐次解；再将激励信号代入方程右端化简得自由项,根据自由项形式求特解；然后根据边界条件求完全解[3].用时域经典法求解例1的基本步骤如下.（1）求齐次解.特征方程为081432=+-αα,可算出41 , 2121==αα.高阶特征根可用MATLAB 的roots 函数计算.齐次解为. 0 , )41()21()(21≥+=n C C n y n n h （2）求方程的特解.将)()43()(n u n x n =代入差分方程右端得自由项为 ⎪⎩⎪⎨⎧≥⋅==-⋅+-1,)43(9130 ,1)1()43(31)()43(1n n n u n u n n n 当1≥n 时,特解可设为n p D n y )43()(=,代入差分方程求得213=D . （3）利用边界条件求完全解.当n =0时迭代求出25)0(=y ,当n ≥1时,完全解的形式为 ,)43(213 )41()21()(21n n n C C n y ⋅++=选择求完全解系数的边界条件可参考文[4]选)1(),0(-y y .根据边界条件求得35,31721=-=C C .注意完全解的表达式只适于特解成立的n 取值范围,其他点要用)(n δ及其延迟表示,如果其值符合表达式则可合并处理.差分方程的完全解为)(])43(213 )41(35)21(317[)1(])43(213 )41(35)21(317[)(25)(n u n u n n y n n n n n n ⋅+⋅+⋅-=-⋅+⋅+⋅-+=δ MATLAB 没有专用的差分方程求解函数,但可调用maple 符号运算工具箱中的rsolve 函数实现[5],格式为y=maple('rsolve({equs, inis},y(n))'),其中:equs 为差分方程表达式, inis 为边界条件,y(n)为差分方程中的输出函数式.rsolve 的其他格式可通过mhelp rsolve 命令了解.在MATLAB 中用时域经典法求解例1中的全响应和单位样值响应的程序如下.clc;clear;format compact;yn=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=(3/4)^n+1/3*(3/4)^(n-1),y(0)=5/2,y(-1)=4},y(n))'),hn=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=0,y(0)=1,y(1)=13/12},y(n))'),3 双零法根据双零响应的定义,按时域经典法的求解步骤可分别求出零输入响应和零状态响应.理解了双零法的求解原理和步骤,实际计算可调用rsolve 函数实现.yzi=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=0,y(-1)=4, y(-2)=12},y(n))'),yzs=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=(3/4)^n+1/3*(3/4)^(n-1),y(0)=1,y(-1)=0},y(n))'),4 变换域法设差分方程的一般形式为)()(00r n x b k n y a r Mr k N k -=-∑∑==.对差分方程两边取单边z 变换,并利用z 变换的位移公式得])()([])()([1010m r m r r M r l k l k k N k z m x z X z b z l y z Y z a ---=-=---=-=∑∑∑∑+=+整理成)()()()()()(00z X z X z B z Y z Y z A +=+形式有. )(, )(110110M M N N z b z b b z B z a z a a z A ----+++=+++=. )()(, )()(110110∑∑∑∑=--=--=--=--==M r r m m r r N k k l l k k z m x b s X zl y a s Y可以看出,由差分方程可直接写出 )(z A 和 )(z B ,系统函数)(/)()(z A z B z H =,将系统函数进行逆z 变换可得单位样值响应.由差分方程的初始状态可算出 )(0z Y ,由激励信号的初始状态可算出 )(0z X ,将激励信号进行z 变换可得 )(z X ,求解z 域代数方程可得输出信号的象函数 , )()()()()()(00z A z Y z X z X z B z Y -+= 对输出象函数进行逆z 变换可得输出信号的原函数)(n y .利用z 变换求解差分方程各响应的步骤可归纳如下:（1）根据差分方程直接写出 )(z A 、 )(z B 和)(z H ,)(z H 的逆变换即为单位样值响应；（2）根据激励信号算出 )(z X ,如激励不是因果序列则还要算出前M 个初始状态值；（3）根据差分方程的初始状态 )(, ),2( ),1(N y y y -⋅⋅⋅--和激励信号的初始状态 )(, ),2( ),1(M x x x -⋅⋅⋅--算出 )(0z Y 和 )(0z X ；（4）在z 域求解代数方程)()()()()()(00z X z X z B z Y z Y z A +=+得输出象函数 )(z Y , )(z Y 的逆变换即为全响应；（5）分析响应象函数的极点来源及在z 平面中的位置,确定自由响应与强迫响应,或瞬态响应与稳态响应；（6）根据零输入响应和零状态响应的定义,在z 域求解双零响应的象函数,对双零响应的象函数进行逆z 变换,得零输入响应和零状态响应.用变换域法求解例1的基本过程如下. 根据差分方程直接写出2181431 )(--+-=z z z A ,1311 )(-+=z z B .系统函数的极点为41,21. 对激励信号进行z 变换得)43/( )(-=z z z X .激励象函数的极点为3/4. 根据差分方程的初始状态算出102123 )(-+-=z z Y .根据激励信号的初始状态算出 0)(0=z X . 对z 域代数方程求解,得全响应的象函数)323161123/()83243125( )(2323-+-+-=z z z z z z z Y . 进行逆z 变换得全响应为)(])43(213 )41(35)21(317[)(n u n y n n n ⋅+⋅+⋅-= 其中,与系统函数的极点对应的是自由响应;与激励象函数的极点对应的是强迫响应. )(z Y 的极点都在z 平面的单位圆内故都是瞬态响应.零输入响应和零状态响应可按定义参照求解.上述求解过程可借助MATLAB 的符号运算编程实现.实现变换域法求解差分方程的m 程序如下: clc;clear;format compact;syms z n %定义符号对象% 输入差分方程、初始状态和激励信号%a=[1,-3/4,1/8],b=[1,1/3], %输入差分方程系数向量y0=[4,12],x0=[0], %输入初始状态,长度分别比a 、b 短1,长度为0时用[]xn=(3/4)^n, %输入激励信号,自动单边处理,u(n)可用1^n 表示% 下面是变换域法求解差分方程的通用程序，极点为有理数时有解析式输出 %N=length(a)-1;M=length(b)-1;%计算长度Az=poly2sym(a,'z')/z^N;Bz=poly2sym(b,'z')/z^M;%计算A(z)和B(z)Hz=Bz/Az;disp('系统函数H(z):'),sys=filt(b,a),%计算并显示系统函数hn=iztrans(Hz);disp('单位样值响应h(n)='),pretty(hn),%计算并显示单位样值响应Hzp=roots(a);disp('系统极点:');Hzp,%计算并显示系统极点Xz=ztrans(xn);disp('激励象函数X(z)='),pretty(Xz),%激励信号的单边z 变换Y0z=0;%初始化Y0(z),求Y0(z)注意系数标号与变量下标的关系for k=1:N;for l=-k:-1;Y0z = Y0z+a(k+1)*y0(-l)*z^(-k-l);endenddisp('初始Y0(z)'),Y0z,%系统初始状态的z 变换X0z=0;%初始化X0(z),求X0(z)注意系数标号与变量下标的关系for r=1:M;for m=-r:-1;X0z = X0z+b(r+1)*x0(-m)*z^(-r-m);endenddisp('初始X0(z)'),X0z,%激励信号起始状态的z 变换Yz=(Bz*Xz+X0z-Y0z)/Az;disp('全响应的z 变换Y(z)'),pretty(simple(Yz)),yn=iztrans(Yz);disp('全响应y(n)='),pretty(yn),% 计算并显示全响应Yziz=-Y0z/Az;disp('零输入象函数Yzi(z)='),pretty(Yziz),%零激励响应的z 变换yzin=iztrans(Yziz);disp('零输入响应yzi(n)='),pretty(yzin),% 计算并显示零输入响应 Yzsz=(Bz*Xz+X0z)/Az;disp('零状态象函数Yzs(z)='),pretty(Yzsz),%零状态响应的z 变换yzsn=iztrans(Yzsz);disp('零状态响应yzs(n)='),pretty(yzsn),% 计算并显示零状态响应该程序的运行过程与手算过程对应,显示在命令窗的运行结果与手算结果相同.。

差分方程的解法分析及MATLAB实现

差分方程的解法分析及MATLAB实现差分方程是描述离散时序系统行为的数学工具。

在离散时间点上，系统的行为由差分方程给出，这是一个递归方程，其中当前时间点的状态取决于之前的状态和其他外部因素。

解差分方程的方法可以分为两类：直接解法和转化为代数方程的解法。

直接解法通过求解差分方程的递归形式来得到解析或数值解。

转化为代数方程的解法则将差分方程转化为代数方程进行求解。

一、直接解法的步骤如下：1.将差分方程表示为递归形式，即将当前时间点的状态表示为之前时间点的状态和其他外部因素的函数。

2.根据初始条件，确定初始时间点的状态。

3.根据递归形式，计算出后续时间点的状态。

以下是一个简单的差分方程的例子：y(n)=2y(n-1)+1，其中n为时间点。

按照上述步骤求解该差分方程：1.将差分方程表示为递归形式：y(n)=2y(n-1)+12.根据初始条件，假设y(0)=1，确定初始时间点的状态。

3.根据递归形式，计算出后续时间点的状态：y(1)=2y(0)+1=2*1+1=3y(2)=2y(1)+1=2*3+1=7y(3)=2y(2)+1=2*7+1=15...依此类推计算出所有时间点的状态。

二、转化为代数方程的解法的步骤如下：1.假设差分方程的解具有指数形式，即y=r^n，其中r为待定参数。

2.将差分方程代入上述假设中，得到r的方程。

3.解得r的值后，再根据初始条件求解出常数值。

4.得到差分方程的解析解。

以下是一个复杂一些的差分方程的例子：y(n)=2y(n-1)+3y(n-2)，其中y(0)=1，y(1)=2按照上述步骤求解该差分方程：1.假设差分方程的解具有指数形式：y=r^n。

2.代入差分方程得到：r^n=2r^(n-1)+3r^(n-2)。

3.整理得到：r^2-2r-3=0。

4.解得r的值为：r1=-1，r2=35.根据初始条件求解出常数值：y(0)=c1+c2=1，y(1)=c1-c2=2、解得c1=1.5，c2=-0.56.得到差分方程的解析解：y(n)=1.5*(-1)^n+-0.5*3^n。

怎样用Matlab求解差分方程题解读

利用plot 绘图观察数量变化趋势
可以用不同线型和颜色绘图 r g b c m y k w 分别表示红绿兰兰绿洋红黄黑白色： + o * . X s d 表示不同的线型
plot(k,y1,k,y2,k,y3) 在同一坐标系下画图 plot(k,y2,':')
>> plot(k,y2,'--')
Matlab实现
首先建立一个关于变量n ，r的函数 function x=sqh(n,r) a=1+r; x=100; for k=1:n x(k+1)=a*x(k); end
在command窗口里调用sqh函数
k=(0:20)';
>> y1=sqh(20,0.0194); >> y2=sqh(20,-0.0324); >> y3=sqh(20,-0.0382); >> round([k,y1',y2',y3'])
将种群按年龄等间隔的分成若干个年龄组，时间也离散化为时段，给定各年龄组种群的繁殖率和死亡率，建立按年龄分组的种群增长模型，预测未来各年龄组的种群数量，并讨论时间充分长以后的变化趋势。
模型及其求解设种群按年龄等间隔的分成n个年龄组，记 i=1,2,· · · ，n,时段记作k=0,1,2· · · ,且年龄组区间与时段长度相等(若5岁为一个年龄组，则5年为一个时段)。以雌性个体为研究对象记在时段k第i年龄组的数量为xi(k);第i年龄组的繁殖率为bi，表示每个个体在一个时段内繁殖的数量；第i年龄组死亡率为di，表示一个时段内死亡数与总数的比，si=1-di是存活率。
1,2 1, xk 0( k )

用Matlab求解差分方程问题

高阶线性常系数差分方程
如果第k+1时段变量Xk+1不仅取决于第k时段变量Xk，而且与以前时段变量有关，就要用高阶差分方程来描述
一年生植物的繁殖
一年生植物春季发芽，夏天开花，秋季产种，没有腐烂，风干，被人为掠取的那些种子可以活过冬天，其中一部分能在第2年春季发芽，然后开花，产种，其中的另一部分虽未能发芽，但如又能活过一个冬天，则其中一部分可在第三年春季发芽，然后开花，产种，如此继续，一年生植物只能活1年，而近似的认为，种子最多可以活过两个冬天，试建立数学模型研究这种植物数量变化的规律，及它能一直繁殖下去的条件。
1, 2

5 10 2
b
植物能一直繁殖下去的条件是b>0.191
线性常系数差分方程组
汽车租赁公司的运营
一家汽车租赁公司在3个相邻的城市运营，为方便顾客起见公司承诺，在一个城市租赁的汽车可以在任意一个城市归还。根据经
验估计和市场调查，一个租赁期内在A市租赁的汽车在A,B,C市归还的比例分别为0.6,0.3,0.1;在B市租赁的汽车归还比例0.2,0.7,0.1;C市租赁的归还比例分别为0.1,0.3,0.6。若公司开业时将600辆汽车平均分配到3个城市，建立运营过程中汽车数量在3个城市间转移的模型，并讨论时间充分长以后的变化趋势。

x2
(k
1)

0.3
0.7
0.3

x2
(
k
)

x3 (k 1) 0.1 0.1 0.6 x3 (k )
function x=czqc(n) A=[0.6,0.2,0.1;0.3,0.7,0.3;0.1,0.1,0.6]; x(:,1)=[200,200,200]'; for k=1:n

MATLAB中的差分方程建模与求解方法

MATLAB中的差分方程建模与求解方法引言差分方程是数学中常见的一种方程类型，是一种离散形式的微分方程。

在实际问题中，差分方程能够提供对系统的离散描述，对于动态模型的建立和求解具有重要作用。

MATLAB作为一种功能强大的数值计算软件，其内置了丰富的工具箱和函数，为差分方程的建模和求解提供了便利。

一、差分方程的建模差分方程的建模是将实际问题转化为数学方程的过程。

在MATLAB中，差分方程的建模可以通过定义离散系统的状态和状态转移方程来实现。

下面以一个简单的例子说明差分方程的建模过程。

假设有一个人口增长模型，人口数在每年增加10%，则差分方程可以表示为：P(n+1) = P(n) + 0.1 * P(n)，其中P(n)表示第n年的人口数，P(n+1)表示第n+1年的人口数。

在MATLAB中，可以通过定义一个函数来描述差分方程的状态转移方程，代码如下：```matlabfunction Pn = population_growth(Pn_minus_1)growth_rate = 0.1;Pn = Pn_minus_1 + growth_rate * Pn_minus_1;end```上述代码定义了一个名为"population_growth"的函数，该函数的输入参数为上一年的人口数"Pn_minus_1"，输出为当前年的人口数"Pn"。

其中，growth_rate表示人口增长率，根据差分方程的定义，将上一年的人口数乘以增长率再加上本身，即可得到当前年的人口数。

二、差分方程的求解方法在MATLAB中，差分方程的求解可以通过多种方法实现。

下面介绍两种常用的差分方程求解方法：欧拉法和四阶龙格-库塔法。

1. 欧拉法（Euler's method）欧拉法是差分方程求解中最简单直观的一种方法。

其基本思想是通过离散化的方式逐步逼近连续函数的解。

具体步骤如下：1) 将时间段分割成若干个小区间；2) 根据差分方程的状态转移方程，在每个小区间内进行计算；3) 迭代计算直到达到指定的时间点。

matlab解差分方程

一阶线性常系数差分方程的解、一阶线性常系数差分方程的解、平衡点及其稳定性
x
k +1
= ax
k
+ b
• 自然环境下,b=0 • 人工孵化条件下
x
k
= a
k
x
0
x k = a k x0 + b (1 + a + ⋯ + a k −1 )
= a
k
x
0
1 − a k + b 1 − a
1 b − a
• 令xk=xk+1=x得
Matlab实现
• • • • • • • 首先建立一个关于变量n ，r的函数 function x=sqh(n,r) a=1+r; x=100; for k=1:n x(k+1)=a*x(k); end
• 在command窗口里调用sqh函数 k=(0:20)'; >> y1=sqh(20,0.0194); >> y2=sqh(20,-0.0324); >> y3=sqh(20,-0.0382); >> round([k,y1',y2',y3'])
• 以k=0时x0=M代入，递推n次可得n年后本息为
xn = (1 + r ) M
n
• 污水处理厂每天可将处理池的污水浓度降低一个固定比例q，问多长时间才能将污水浓度降低一半？ • 记第k天的污水浓度为ck,则第k+1天的污水浓度为 ck+1=(1-q)ck,k=0,1,2,···· 从k=0开始递推n次得
x1 ( k + 1) = 0.6 x1 ( k ) + 0.2 x2 ( k ) + 0.1x3 ( k ) x2 ( k + 1) = 0.3 x1 ( k ) + 0.7 x2 ( k ) + 0.3 x3 ( k ) x ( k + 1) = 0.1x ( k ) + 0.1x ( k ) + 0.6 x ( k ) 1 2 3 3λ1 2=5 ± 2

用Matlab求解差分方程问题

Xk-1决定的部分是 a1bcXk-1,由Xk-2决定的部分是 a2b(1-a1)bcXk-2
Xk= a1bcXk-1 + a2b(1-a1)bcXk-2
Xk= a1bcXk-1 + a2b(1-a1)bcXk-2
• 实际上，就是Xk= pXk-1 + qXk-2 我们需要知道x0,a1,a2,c, 考察b不同时，种子繁殖的情况。在这里假设 • X0=100,a1=0.5,a2=0.25,c=10,b=0.18~0.20 • 这样可以用matlab计算了
Xk= a1bcXk-1 + a2b(1-a1)bcXk-2
• • • • • • • • Function x=zwfz(x0,n,b) C=10;a1=0.5;a2=0.25; p=a1*b*c;q=a2*b*(1-a1)*b*c; X1=x0; X2=p*(x1); for k=3:n X(k)=p*(xk-1)+q*(xk-2); end
k=0:10; plot(k,x) ，grid gtext('x1(k)'), gtext('x2(k)'), gtext('x3(k)')
280 260 240 220 200 180 160 140 120
x2(k)
x1(k)
x3(k) 0 1 2 3 4 5 6 7 8 9 10
• 可以看到时间充分长以后3个城市汽车数量趋于180，300，120 • 可以考察这个结果与初始条件是否有关 • 若最开始600辆汽车都在A市，可以看到变化时间充分长以后，各城市汽车数量趋于稳定，与初始值无关
直接输入x（:,1）的值即可
x(:,1)=[600,0,0]; round(x'); plot(k,x),grid

matlab差分法解偏微分方程

Matlab 差分法解偏微分方程1.引言解偏微分方程是数学和工程领域中的一项重要课题，它在科学研究和工程实践中具有广泛的应用。

而 Matlab 差分法是一种常用的数值方法，用于求解偏微分方程。

本文将介绍 Matlab 差分法在解偏微分方程中的应用，包括原理、步骤和实例。

2. Matlab 差分法原理差分法是一种离散化求解微分方程的方法，通过近似替代微分项来求解微分方程的数值解。

在 Matlab 中，差分法可以通过有限差分法或者差分格式来实现。

有限差分法将微分方程中的导数用有限差分替代，而差分格式指的是使用不同的差分格式来近似微分方程中的各个项，通常包括前向差分、后向差分和中心差分等。

3. Matlab 差分法步骤使用 Matlab 差分法解偏微分方程一般包括以下步骤：（1）建立离散化的区域：将求解区域离散化为网格点或节点，并确定网格间距。

（2）建立离散化的时间步长：对于时间相关的偏微分方程，需要建立离散化的时间步长。

（3）建立离散化的微分方程：使用差分法将偏微分方程中的微分项转化为离散形式。

（4）建立迭代方程：根据离散化的微分方程建立迭代方程，求解数值解。

（5）编写 Matlab 代码：根据建立的迭代方程编写 Matlab 代码求解数值解。

（6）求解并分析结果：使用 Matlab 对建立的代码进行求解，并对结果进行分析和后处理。

4. Matlab 差分法解偏微分方程实例假设我们要使用 Matlab 差分法解决以下一维热传导方程：∂u/∂t = α * ∂^2u/∂x^2其中 u(x, t) 是热传导方程的温度分布，α 是热扩散系数。

4.1. 离散化区域和时间步长我们将求解区域离散化为网格点，分别为 x_i，i=1,2,...,N。

时间步长为Δt。

4.2. 离散化的微分方程使用中心差分格式将偏微分方程中的导数项离散化得到：∂u/∂t ≈ (u_i(t+Δt) - u_i(t))/Δt∂^2u/∂x^2 ≈ (u_i-1(t) - 2u_i(t) + u_i+1(t))/(Δx)^2代入原偏微分方程可得离散化的微分方程：(u_i(t+Δt) - u_i(t))/Δt = α * (u_i-1(t) - 2u_i(t) + u_i+1(t))/(Δx)^24.3. 建立迭代方程根据离散化的微分方程建立迭代方程：u_i(t+Δt) = u_i(t) + α * Δt * (u_i-1(t) - 2u_i(t) + u_i+1(t))/(Δx)^24.4. 编写 Matlab 代码使用以上建立的迭代方程编写 Matlab 代码求解热传导方程。

matlab解差分方程

解析法
通过数学变换和求解技巧，得到差分方程的解析解。
数值法
通过数值计算方法，如欧拉法、龙格-库塔法等，得到差分方程的近似解。
02
Matlab在解差分方程中的应用
Matlab的符号计算功能
符号计算
Matlab提供了符号计算的功能，可以用于解决差分方程的符号解。通过符号计算，可以找到差分方程的通解或特解，并分析解的性质。
析解。
THANKS
要点二
递推法
通过递推公式求解高阶差分方程，适用于一阶高阶差分方程。
04
差分方程的数值解法
欧拉方法
总结词
简单易行，但精度较低
详细描述
欧拉方法是一种简单的数值方法，用于求解差分方程。它基于差分方程的递推性质，通过迭代的方式逐步逼近解。由于其简单性，欧拉方法在许多情况下是首选的数值方法。然而，由于其精度较低，对于需要高精度解的问题，可能需要采用更高级的数值方法。
特征方程法
通过求解特征方程来求解线性差分方程，适用于多阶线性差分方程。
非线性差分方程的解析解法
迭代法
通过迭代公式求解非线性差分方程，适用于一阶非线性差分方程。
解析法
通过解析表达式求解非线性差分方程，适用于多阶非线性差分方程。
高阶差分方程的解析解法
要点一
降阶法
将高阶差分方程转化为低阶差分方程，再求解低阶差分方程。
05
Matlab实现差分方程的解法示例
解析解法的示例
差分方程
$y(n+1) - y(n) = 0$
解析解
$y(n) = y(0)$
解释
该差分方程表示每一项都与前一项相等，因此解为常数。
数值解法的示例

用matlab解差分方程

• B=1; A=[1,-a]; %差分方程系数
• xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xi
• yn=filter(B,A,xn,xi); %调用filter解差分方程，求系统输出信号y(n)
• n=0:length(yn)-1;
• subplot(3,2,1);stem(n,yn,'.')
• title('(a)'); xlabel('(n)');ylabel('(n)');%ep141.m: 调用filter解差分方程y(n)-ay(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,-a]; %差分方程系数
ห้องสมุดไป่ตู้
• xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xi
• yn=filter(B,A,xn,xi); %调用filter解差分方程，求系统输出信号y(n)
图形表现
用matlab解差分方程
本例子截取书上第19页，
matlab求解程序
• %ep141.m: 调用filter解差分方程y(n)-ay(n-1)=x(n)
• a=0.8; ys=1;
%设差分方程系数a=0.8，初始状态：y(-1)=1
• xn=[1,zeros(1,30)] %x(n)=单位脉冲序列，长度N=31

用matlab求解差分方程

一阶线性常系数差分方程的解、平衡点及其稳定性
01
单击此处添加正文，文字是您思想的提炼，请尽量言简意赅地阐述观点。
自然环境下,b=0
02
单击此处添加正文，文字是您思想的提炼，请尽量言简意赅地阐述观点。
人工孵化条件下
03
差分方程的平衡点
令xk=xk+1=x得
04
单击此处添加正文，文字是您思想的提炼，请尽量言简意赅地阐述观点。
添加标题
X(k)=p*(xk-1)+q*(xk-2);
添加标题
end
K=(0:20)’; Y1=zwfz(100,21,0.18); Y2=zwfz(100,21,0.19); Y3=zwfz(100,21,0,20); Round([k,y1’,y2’,y3’]) Plot(k,y1,k,y2,’:’,k,y3,’o’), Gtext(‘b=0.18’),gtext(‘b=0.19’),gtext(‘b=0.20’)
plot(k,y2,':') >> plot(k,y2,'--') >> plot(k,y2,'r') >> plot(k,y2,'y') >> plot(k,y2,'y',k,y1,':') >> plot(k,y2,k,y1,':') >> plot(k,y2,'oy',k,y1,':') 用gtext(‘r=0.0194’),gtext(‘r=-0.0324’),gtext(‘r=-0.0382’)在图上做标记。
可以看到时间充分长以后3个城市汽车数量趋于180，300，120 可以考察这个结果与初始条件是否有关若最开始600辆汽车都在A市，可以看到变化时间充分长以后，各城市汽车数量趋于稳定，与初始值无关

差分方程的解法分析及其MATLAB实现

差分方程的解法分析及其MATLAB实现张登奇;彭仕玉【摘要】差分方程是描述离散时间系统的数学模型，求解差分方程是分析离散时间系统的重要内容，常用的求解方法有迭代法、时域经典法、双零法和变换域法。

文章根据各种方法的求解原理，分别介绍了不同方法的求解步骤，结合实例列出了这些方法的求解过程及MATLAB实现程序。

%The difference equation is mathematical model to describe discrete-time systems, to solve the differential equation is an important part to analyze discrete-time systems, commonly used methods are: iterative method, classical time-domain method, double zero response method and the transform-domain method. This paper introduces the solving steps of these different methods according to the corresponding principles with examples and lists these corresponding MATLAB programs.【期刊名称】《湖南理工学院学报（自然科学版）》【年(卷),期】2014(000)003【总页数】5页(P28-32)【关键词】离散时间系统;差分方程;时域分析法;变换域分析法;MATLAB【作者】张登奇;彭仕玉【作者单位】湖南理工学院信息与通信工程学院，湖南岳阳 414006;湖南理工学院信息与通信工程学院，湖南岳阳 414006【正文语种】中文【中图分类】TN911.72;O175.7线性常系数差分方程是描述线性时不变离散时间系统的数学模型, 求解差分方程是分析离散时间系统的重要内容. 在《信号与系统》课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域法[1]. 迭代法可手工逐次代入求解, 也可编程用计算机求解, 该方法原理简单, 缺点是只能得到数值解.时域经典法先求齐次解和特解, 再用边界条件确定待定系数得完全解, 该方法数学过程清晰, 但求解过程麻烦. 双零法分别求零输入响应和零状态响应, 再通过叠加得到全响应. 该方法物理意义清晰, 但求解过程依然麻烦. 变换域法利用z变换将差分方程变换成代数方程求解, 该方法简便高效, 是求解差分方程的重要方法. 本文根据不同方法的求解原理, 分别介绍各种方法的求解步骤, 结合实例列出这些方法的求解过程和MATLAB实现程序.差分方程本身就是一个递推方程, 根据初始状态和激励信号依次迭代就可算出输出序列. 迭代法是解差分方程的基础方法, 如果所需输出序列个数较少（如计算边界条件）可手工直算, 如需计算大量输出可利用计算机编程实现.现结合实例介绍迭代法的计算过程.例1已知离散系统的差分方程为激励信号为初始状态为y(−1)=4,y(−2)=12.求系统响应.根据激励信号和初始状态, 手工依次迭代可算出利用MATLAB中的filter函数实现迭代过程的m程序如下:clc;clear;format compact;a=[1, -3/4, 1/8], b=[1, 1/3, 0], %输入差分方程系数向量, 不足补0对齐n=0:10;xn=(3/4).^n, %输入激励信号zx=[0, 0], zy=[4, 12], %输入初始状态zi=filtic(b, a, zy, zx), %计算等效初始条件[yn, zf]=filter(b, a, xn, zi), %迭代计算输出和后段等效初始条件MATLAB提供的filter函数是一个内建函数, 用type命令看不到程序代码.为了理解迭代思想, 下面根据图1所示的直接Ⅰ型结构[2], 重写实现迭代法的m程序. clc;clear;format compact;a=[1, -3/4, 1/8], b=[1, 1/3, 0], %输入差分方程系数向量, 不足补0对齐n=0:10;x=(3/4).^n, %输入激励信号zx=[0, 0], zy=[4, 12], %输入初始状态% 下面是按直接Ⅰ型结构迭代的通用程序 %N=length(a)-1, %计算数据存储长度L=length(x), %计算激励信号长度y=zeros(1, L);%输出信号初始化for i=1:L; %逐个计算输出信号for n=1:N;z(n)=b(n+1)*zx(n)-a(n+1)*zy(n);end %分算过去zz=sum(z);%计算输出中的过去分量y(i)=b(1)*x(i)+zz;% 计算当前输出y(n)for n=N:-1:2, zx(n)=zx(n-1);zy(n)=zy(n-1);end%过去数据下移zx(1)=x(i);zy(1)=y(i);%当前的激励和输出变为过去, 以便算下一个输出end%理解filter函数中zf参数的意义zf=zeros(1, N);%初始化zffor k=1:N; %逐个计算zf参数for n=1:N;z(n)=b(n+1)*zx(n)-a(n+1)*zy(n);end %算z(n)zf(k)=sum(z);% 计算第k个zffor k=N:-1:2, zx(k)=zx(k-1);zy(k)=zy(k-1);end;%过去数据下移zx(1)=0;zy(1)=0;% 没有当前的激励和输出变为过去, 算下一个zfendy, zf, %显示输出和zf参数该程序采用直接Ⅰ型结构实现迭代计算, 与filter函数计算的结果相同, zf参数可用于激励信号分段处理时，调用filter函数计算下一段输出所需的等效初始条件.用时域经典法求解差分方程与高等数学中求解微分方程的过程类似: 先求齐次解; 再将激励信号代入方程右端化简得自由项, 根据自由项形式求特解; 然后根据边界条件求完全解[3]. 用时域经典法求解例1的基本步骤如下.(1) 求齐次解. 特征方程为可算出高阶特征根可用MATLAB的roots函数计算. 齐次解为(2) 求方程的特解. 将代入差分方程右端得自由项为当n≥1时, 特解可设为代入差分方程求得(3) 利用边界条件求完全解.当n=0时迭代求出当n≥1时, 完全解的形式为选择求完全解系数的边界条件可参考文[4]选y(0),y(−1). 根据边界条件求得注意完全解的表达式只适于特解成立的n取值范围, 其他点要用δ(n)及其延迟表示, 如果其值符合表达式则可合并处理.差分方程的完全解为MATLAB没有专用的差分方程求解函数, 但可调用maple符号运算工具箱中的rsolve函数实现[5], 格式为y=maple('rsolve({equs, inis}, y(n))'), 其中equs为差分方程表达式, inis为边界条件,y(n)为差分方程中的输出函数式. rsolve的其他格式可通过mhelp rsolve命令了解. 在MATLAB中用时域经典法求解例1中的全响应和单位样值响应的程序如下:clc;clear;format compact;yn=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=(3/4)^n+1/3*(3/4)^(n-1),y(0)=5/2, y(-1)=4}, y(n))'),hn=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=0, y(0)=1, y(1)=13/12},y(n))'),如同时域经典法一样, 应用rsolve函数时要特别注意边界条件的选择问题，且边界条件要连续取点.本例求单位样值响应时的边界条件要选用y(0)和y(1)[6], 在命令窗中显示的结果也要分析n的取值范围.双零法是将完全解分解成物理意义明显的零输入响应和零状态响应分别计算. 零输入响应是激励为零, 由系统的初始状态所产生的响应; 零输入响应要求差分方程右端为零, 故特解为零; 完全解为齐次解形式, 系数可直接由初始状态确定. 零状态响应是初始状态为零, 由激励信号所产生的响应. 计算零状态响应可用时域经典法, 也可用卷积法. 根据双零响应的定义, 按时域经典法的求解步骤可分别求出零输入响应和零状态响应. 理解了双零法的求解原理和步骤, 实际计算可调用rsolve函数实现:yzi=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=0, y(-1)=4, y(-2)=12}, y(n))'), yzs=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=(3/4)^n+1/3*(3/4)^(n-1), y(0)=1, y(-1)=0}, y(n))'),边界条件和n的取值范围这里不再赘述, 各响应的计算结果在命令窗都有显示, 零输入响应和零状态响应之和与前面计算的全响应相等.变换域法是以z变换为数学工具, 先将时域差分方程变换成z域代数方程, 再在z 域求解响应的象函数, 最后将响应象函数逆变换成时域原函数, 是一种间接求解差分方程的计算方法.z变换的线性和位移性是变换域法求解差分方程的重要性质, 这里只列出双边序列的单边z变换位移公式[1]:设差分方程的一般形式为对差分方程两边取单边z变换, 并利用z变换的位移公式, 得整理成A(z)Y(z)+Y0(z)=B(z)X(z)+X0(z)形式, 有可以看出, 由差分方程可直接写出A(z)和B(z), 系统函数将系统函数进行逆z变换可得单位样值响应. 由差分方程的初始状态可算出Y0(z), 由激励信号的初始状态可算出X0(z), 将激励信号进行z变换可得X(z), 求解z域代数方程可得输出信号的象函数对输出象函数进行逆z变换可得输出信号的原函数y(n). 为了计算简单和便于对结果进行对比分析,例1列举的实例中, 系统函数的极点均为单根有理数极点, 且系统函数H(z)的零极点与激励象函数X(z)的零极点均不相同. 参照s域求解微分方程的方法[7], 利用z变换求解差分方程各响应的步骤可归纳如下:(1) 根据差分方程直接写出A(z),B(z)和H(z),H(z)的逆变换即为单位样值响应;(2) 根据激励信号算出X(z), 如激励不是因果序列则还要算出前M个初始状态值;(3) 根据差分方程的初始状态y(−1),y(−2),⋅⋅⋅,y(−N)和激励信号的初始状态x(−1),x(−2),⋅⋅⋅,x(−M)算出Y0(z)和X0(z);(4) 在z域求解代数方程A(z)Y(z)+Y0(z)=B(z)X(z)+X0(z)得输出象函数Y(z),Y(z)的逆变换即为全响应;(5) 分析响应象函数的极点来源及在z平面中的位置, 确定自由响应与强迫响应, 或瞬态响应与稳态响应;(6) 根据零输入响应和零状态响应的定义, 在z域求解双零响应的象函数, 对双零响应的象函数进行逆z变换, 得零输入响应和零状态响应.用变换域法求解例1的基本过程如下.根据差分方程直接写出系统函数的极点为对激励信号进行z变换得激励象函数的极点为根据差分方程的初始状态算出根据激励信号的初始状态算出X0(z)=0.对z域代数方程求解, 得全响应的象函数进行逆z变换得全响应为其中, 与系统函数的极点对应的是自由响应; 与激励象函数的极点对应的是强迫响应.Y(z)的极点都在z平面的单位圆内故都是瞬态响应. 零输入响应和零状态响应可按定义参照求解.上述求解过程可借助MATLAB的符号运算编程实现. 实现变换域法求解差分方程的m程序如下:clc;clear;format compact;syms z n %定义符号对象% 输入差分方程、初始状态和激励信号%a=[1, -3/4, 1/8], b=[1, 1/3], %输入差分方程系数向量y0=[4, 12], x0=[0], %输入初始状态, 长度分别比a、b短1, 长度为0时用[]xn=(3/4)^n, %输入激励信号, 自动单边处理, u(n)可用1^n表示% 下面是变换域法求解差分方程的通用程序，极点为有理数时有解析式输出 % N=length(a)-1;M=length(b)-1;%计算长度Az=poly2sym(a, 'z')/z^N;Bz=poly2sym(b, 'z')/z^M;%计算A(z)和B(z)Hz=Bz/Az;disp('系统函数H(z):'), sys=filt(b, a), %计算并显示系统函数hn=iztrans(Hz);disp('单位样值响应h(n)='), pretty(hn), %计算并显示单位样值响应Hzp=roots(a);disp('系统极点:');Hzp, %计算并显示系统极点Xz=ztrans(xn);disp('激励象函数X(z)='), pretty(Xz), %激励信号的单边z变换Y0z=0;%初始化Y0(z), 求Y0(z)注意系数标号与变量下标的关系for k=1:N;for l=-k:-1;Y0z = Y0z+a(k+1)*y0(-l)*z^(-k-l);endenddisp('初始Y0(z)'), Y0z, %系统初始状态的z变换X0z=0;%初始化X0(z), 求X0(z)注意系数标号与变量下标的关系for r=1:M;for m=-r:-1;X0z = X0z+b(r+1)*x0(-m)*z^(-r-m);endenddisp('初始X0(z)'), X0z, %激励信号起始状态的z变换Yz=(Bz*Xz+X0z-Y0z)/Az;disp('全响应的z变换Y(z)'), pretty(simple(Yz)),yn=iztrans(Yz);disp('全响应y(n)='), pretty(yn), % 计算并显示全响应Yziz=-Y0z/Az;disp('零输入象函数Yzi(z)='), pretty(Yziz), %零激励响应的z变换yzin=iztrans(Yziz);disp('零输入响应yzi(n)='), pretty(yzin), % 计算并显示零输入响应Yzsz=(Bz*Xz+X0z)/Az;disp('零状态象函数Yzs(z)='), pretty(Yzsz), %零状态响应的z变换yzsn=iztrans(Yzsz);disp('零状态响应yzs(n)='), pretty(yzsn), % 计算并显示零状态响应该程序的运行过程与手算过程对应, 显示在命令窗的运行结果与手算结果相同.求解差分方程有多种方法. 迭代法计算原理简单, 调用filter函数时要注意初始状态的等效变换. 经典法数学思路清晰, 调用rsolve函数时要注意边界条件的选择问题. 双零法物理意义清楚, 可根据定义用经典法计算. 变换域法简便高效, 适合计算各类响应. 本文列举的变换域法求解实例可作为课程教学的补充材料, 编写的变换域法求解程序对教学科研也有一定的应用价值.【相关文献】[1] 郑君里, 应启珩, 杨为理. 信号与系统[M]. 第2版. 北京: 高等教育出版社, 2000:16, 64[2] 高西全, 丁玉美. 数字信号处理[M]. 第3版. 西安: 西安电子科技大学出版社, 2008:129[3] 陈从颜, 翟军勇. 信号与系统基础[M]. 北京: 机械工业出版社, 2009:19[4] 朱玲赞. 差分方程的边界条件和离散系统的初始状态[J]. 电气电子教学学报, 2001(03): 35～37[5] 张志勇. 精通MATLAB6.5[M]. 北京: 北京航空航天大学出版社, 2003:229[6] 史林杰. 离散时间系统边界条件的确定准则[J]. 电工教学, 1997(02): 22～24[7] 张登奇, 张璇. 连续时间系统的s域分析及MATLAB实现[J]. 湖南理工学院学报(自然科学版), 2012(02): 26～29。

matlab有限差分法

matlab有限差分法一、前言Matlab是一种广泛应用于科学计算和工程领域的计算机软件，它具有简单易学、功能强大、易于编程等优点。

有限差分法（Finite Difference Method）是一种常用的数值解法，它将微分方程转化为差分方程，通过对差分方程进行离散化求解，得到微分方程的数值解。

本文将介绍如何使用Matlab实现有限差分法。

二、有限差分法基础1. 有限差分法原理有限差分法是一种通过将微分方程转化为离散形式来求解微分方程的数值方法。

其基本思想是将求解区域进行网格划分，然后在每个网格点上进行逼近。

假设要求解一个二阶常微分方程：$$y''(x)=f(x,y(x),y'(x))$$则可以将其转化为离散形式：$$\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}=f(x_i,y_i,y'_i)$$其中$h$为网格步长，$y_i$表示在$x_i$处的函数值。

2. 一维情况下的有限差分法对于一维情况下的常微分方程：$$\frac{d^2 y}{dx^2}=f(x,y,y')$$可以使用中心差分法进行离散化：$$\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}=f(x_i,y_i,y'_i)$$这个方程可以写成矩阵形式：$$A\vec{y}=\vec{b}$$其中$A$为系数矩阵，$\vec{y}$为函数值向量，$\vec{b}$为右端项向量。

三、Matlab实现有限差分法1. 一维情况下的有限差分法假设要求解的方程为：$$\frac{d^2 y}{dx^2}=-\sin(x)$$首先需要确定求解区域和网格步长。

在本例中，我们将求解区域设为$[0,2\pi]$，网格步长$h=0.01$。

则可以通过以下代码生成网格：```matlabx = 0:0.01:2*pi;```接下来需要构造系数矩阵和右端项向量。

根据上面的公式，系数矩阵应该是一个三对角矩阵，可以通过以下代码生成：```matlabn = length(x)-2;A = spdiags([-ones(n,1), 2*ones(n,1), -ones(n,1)], [-1 0 1], n, n); ```其中`spdiags`函数用于生成一个稀疏矩阵。

一维流体差分 matlab

一维流体差分 matlab
在处理一维流体动力学问题时，差分方法是一种常用的数值求
解技术。

在Matlab中，你可以使用差分方法来离散化一维流体动力
学方程，然后求解离散化后的方程以获得数值解。

首先，你需要将一维流体动力学方程离散化。

例如，如果你正
在处理一维不可压缩流体的流动问题，你可以使用连续方程和动量
方程来描述流体的行为。

然后，你可以将这些偏微分方程转化为差
分方程，例如使用有限差分方法。

在Matlab中，你可以编写一个函数来表示差分方程，并使用内
置的数值求解器（如ode45）来求解这些方程。

你需要定义网格、
边界条件和初始条件，并使用适当的差分格式（如向前差分、向后
差分或中心差分）来离散化偏微分方程。

然后，你可以使用Matlab
的矩阵运算和求解器来解决离散化后的方程。

此外，Matlab还提供了一些专门用于求解偏微分方程的工具箱，如Partial Differential Equation Toolbox，它包含了一些内置
的函数和工具，可以帮助你更轻松地处理一维流体动力学问题的数
值求解。

总之，使用Matlab进行一维流体动力学问题的差分求解涉及到离散化方程、定义边界条件和初始条件、选择合适的差分格式，以及使用Matlab的数值求解器或工具箱来求解离散化后的方程。

希望这些信息能够帮助你更好地理解在Matlab中应用差分方法求解一维流体动力学问题的过程。