下载文档原格式
第8章 常微分方程的数值求解所谓的常微分方程就是把自变量t 和它的函数y 以及它的微商dy/dt 、d 2y/dt 2、…d n y/dt n 相联系的一个关系式0,...,,,22=⎪⎪⎭⎫ ⎝⎛n n dt y d dt y d dt dy y t f (1) 一个微分方程不只有一个或几个解，而是有无数个一簇解。
如一阶常微分方程dy/dt=1-e -t 的解为y(t)=t+e -t +C 。
C 为积分常数，C 取任意数值时，函数y(t)都满足微分方程。
因此，解有无数个。
如图在实际应用中，并不要求把所有的解都求出来，而是求满足某种指定条件的解。
这个条件通常称定解条件。
一个最重要的定解条件是初值条件。
对于上述方程，初值条件是：()00y t y =，()'00y dt t dy =，())2(0202y dt t dy =，…，())1(0101---=n n n y dt t dy (2) x 0是自变量的某个指定的“初值”，而0y 、'0y 、)2(0y 、…、)1(0-n y 则是未知函数及其到n －1阶微商的指定“初值”。
求解满足这样初值条件的微分方程问题称为初值问题。
如果能从方程(1)将n n dt y d 解出，则微分方程变成下面的形式⎪⎪⎭⎫ ⎝⎛=--1122,...,,,n n n n dt y d dt y d dt dy y t f dt y d (3) 这里的f 与式(1)中的f 不同，它是n+1个自变量的已知函数。
这种微分方程称为正规形微分方程。
而式(1)有时称为隐微分方程。
我们只考虑正规形微分方程，而且是一阶常微分方程的初值问题： ()y t f dtdy,=，()00y t y =， (4) 设函数()y t f ,在区域T t t ≤≤0，∞<u 内连续，并且存在常数L(Lipschitz 常数)，对所有],[00T t t ∈和y 1、y 2，有()()2121,,y y L y t f y t f -≤-，由常微分方程理论得知，初值问题(4)在区间[t 0, T]有唯一解，且连续可微。
MATLAB是一种用于科学计算和工程应用的高级编程语言和交互式环境。
它在数学建模、模拟和分析等方面有着广泛的应用。
在MATLAB 中，常微分方程的数值求解是一个常见的应用场景。
在实际工程问题中，通常需要对常微分方程进行数值求解来模拟系统的动态行为。
本文将介绍MATLAB中对常微分方程进行数值求解的快速方法。
1. 基本概念在MATLAB中，可以使用ode45函数来对常微分方程进行数值求解。
ode45是一种常用的Runge-Kutta法，它可以自适应地选取步长，并且具有较高的数值精度。
使用ode45函数可以方便地对各种类型的常微分方程进行求解，包括一阶、高阶、常系数和变系数的微分方程。
2. 函数调用要使用ode45函数进行常微分方程的数值求解，需要按照以下格式进行函数调用：[t, y] = ode45(odefun, tspan, y0)其中，odefun表示用于描述微分方程的函数，tspan表示求解的时间跨度，y0表示初值条件，t和y分别表示求解得到的时间序列和对应的解向量。
3. 示例演示为了更好地理解如何使用ode45函数进行常微分方程的数值求解，下面我们以一个具体的例子来进行演示。
考虑如下的一阶常微分方程：dy/dt = -2*y其中，y(0) = 1。
我们可以编写一个描述微分方程的函数odefun：function dydt = odefun(t, y)dydt = -2*y;按照上述的函数调用格式，使用ode45函数进行求解：tspan = [0 10];y0 = 1;[t, y] = ode45(odefun, tspan, y0);绘制出解曲线：plot(t, y);4. 高级用法除了基本的函数调用方式外，MATLAB中还提供了更多高级的方法来对常微分方程进行数值求解。
可以通过设定选项参数来控制数值求解的精度和稳定性，并且还可以对刚性微分方程进行求解。
5. 性能优化在实际工程应用中，常常需要对大规模的常微分方程进行数值求解。
第8章 MATLAB方程数值求解习题8一、选择题1．下列方法中与线性方程组求解无关的是（）。
CA．左除B．矩阵求逆C．矩阵转置D．矩阵分解2．对于系数矩阵A的阶数很大，且零元素较多的大型稀疏矩阵线性方程组，非常适合采用（）求解。
BA．直接法B．迭代法C．矩阵求逆D．左除3．已知函数文件fx.m：function f=fx(x)f=2*x.^2+5*x-1;则求f(x)=2x2+5x-1=0在x0=-2附近根的命令是（）。
DA．z=fzero(fx,0.5) B．z=fzero(@fx,0.5)C．z=fzero(fx,-2); D．z=fzero(@fx,-2);4．已知：fx=@(x) 2*x.^2+5*x-1;则求f(x)=2x2+5x-1=0在x0=-2附近根的命令是（）。
CA．z=fzero(fx,0.5) B．z=fzero(@fx,0.5)C．z=fzero(fx,-2); D．z=fzero(@fx,-2);5．下列选项中不能用于求常微分方程数值解的函数是（）。
AA．ode10 B．ode23 C．ode45 D．ode113二、填空题1．线性方程组的求解方法可以分为两类，一类是，另一类是。
前者是在没有舍入误差的情况下，通过有限步的初等运算来求得方程组的解；后者是先给定一个解的，然后按照一定的算法不断用变量的旧值递推出新的值。
直接法，迭代法，初始值2．MA TLAB用函数来求单变量非线性方程的根。
对于非线性方程组，则用函数求其数值解。
fzero，fsolve3．用数值方法求解常微分方程的初值问题，一般都是用系列函数，包括ode23、ode45等函数，各有不同的适用场合。
ode4．ode23、ode45等函数是针对一阶常微分方程组的，对于高阶常微分方程，需先将它转化为一阶常微分方程组，即。
状态方程三、应用题21．分别用矩阵除法以及矩阵分解求线性方程组的解。
⎪⎩⎪⎨⎧=+-=++=++57347310532)1(z y x z y x z y x 123134124345132(2)53241x x x x x x x x x x x +-=⎧⎪+-=⎪⎨--+=⎪⎪+=-⎩(1): 矩阵除法:A=[2,3,5;3,7,4;1,-7,1]; B=[10,3,5];%B 是行向量 x=A\B'%将B 变成列向量 矩阵分解:A=[2,3,5;3,7,4;1,-7,1]; B=[10,3,5];%B 是行向量 [L,U]=lu(A); x=U\(L\B') (2):和上面的程序一样。
MATLAB是一种流行的数学软件，用于解决各种数学问题，包括微分方程的数值积分。
微分方程是许多科学和工程问题的数学描述方式，通过数值积分可以得到微分方程的数值解。
本文将介绍在MATLAB中如何进行微分方程的数值积分，以及一些相关的技巧和注意事项。
一、MATLAB中微分方程的数值积分的基本方法1. 常微分方程的数值积分在MATLAB中，常微分方程的数值积分可以使用ode45函数来实现。
ode45是一种常用的数值积分函数，它使用4阶和5阶Runge-Kutta 方法来求解常微分方程。
用户只需要将微分方程表示为函数的形式，并且提供初值条件，ode45就可以自动进行数值积分，并得到微分方程的数值解。
2. 偏微分方程的数值积分对于偏微分方程的数值积分，在MATLAB中可以使用pdepe函数来实现。
pdepe可以求解具有定解条件的一维和二维偏微分方程，用户只需要提供偏微分方程的形式和边界条件，pdepe就可以进行数值积分，并得到偏微分方程的数值解。
二、在MATLAB中进行微分方程数值积分的注意事项1. 数值积分的精度和稳定性在进行微分方程的数值积分时，需要注意数值积分的精度和稳定性。
如果数值积分的精度不够，可能会导致数值解的误差过大；如果数值积分的稳定性差，可能会导致数值解发散。
在选择数值积分方法时，需要根据具体的微分方程来选择合适的数值积分方法，以保证数值解的精度和稳定性。
2. 初值条件的选择初值条件对微分方程的数值解有很大的影响，因此在进行微分方程的数值积分时，需要选择合适的初值条件。
通常可以通过对微分方程进行分析，或者通过试验求解来确定合适的初值条件。
3. 数值积分的时间步长在进行微分方程的数值积分时，需要选择合适的时间步长，以保证数值积分的稳定性和效率。
选择时间步长时，可以通过试验求解来确定合适的时间步长，以得到最优的数值解。
三、MATLAB中微分方程数值积分的实例以下通过一个简单的例子来演示在MATLAB中如何进行微分方程的数值积分。
用matlab 求解常微分方程在MATLAB 中，由函数dsolve ()解决常微分方程（组）的求解问题，其具体格式如下：r = dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v')'eq1,eq2,...'为微分方程或微分方程组，'cond1,cond2,...',是初始条件或边界条件，'v'是独立变量，默认的独立变量是't'。
函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。
例1：求解常微分方程1dy dx x y =+的MATLAB 程序为：dsolve('Dy=1/(x+y)','x')，注意，系统缺省的自变量为t ，因此这里要把自变量写明。
其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X 。
例2：求解常微分方程的MATLAB 程序为：2'''0yy y −=Y2=dsolve('y*D2y-Dy^2=0','x')Y2=dsolve('D2y*y-Dy^2=0','x')我们看到有两个解，其中一个是常数0。
例3：求常微分方程组253ttdxx y edtdyx y edt⎧++=⎪⎪⎨⎪−−=⎪⎩通解的MATLAB程序为：[X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t')例4：求常微分方程组2210cos,24,tttdx dyx t xdt dtdx dyy e ydt dt=−=⎧+−==⎪⎪⎨⎪++==⎪⎩2通解的MATLAB程序为：[X,Y]=dsolve('Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2,y(0)=0','t')以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。
第8章MATLAB方程数值求解例8-1用直接解法求解下列线性方程组。
程序如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';x=A\b例8-2用LU分解求解例8-1中的线性方程组。
程序如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[L,U]=lu(A);x=U\(L\b)例8-3 用QR分解求解例8-1中的线性方程组。
程序如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[Q,R]=qr(A);x=R\(Q\b)例8-4 用Cholesky分解求解例8-1中的线性方程组。
命令如下：>> A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; >> b=[13,-9,6,0]';>> R=chol(A)Jacobi迭代法的MA TLAB函数文件jacobi.m如下：function [y,n]=jacobi(A,b,x0,ep)if nargin==3ep=1.0e-6;elseif nargin<3errorreturnendD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=D\(L+U);f=D\b;y=B*x0+f;n=1; %迭代次数while norm(y-x0)>=epx0=y;y=B*x0+f;n=n+1;end例8-5 用Jacobi迭代法求解下列线性方程组。
设迭代初值为0，迭代精度为10-6。
在程序中调用函数文件jacobi.m，程序如下：A=[10,-1,0;-1,10,-2;0,-2,10];b=[9,7,6]';[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)Gauss-Serdel迭代法的MA TLAB函数文件gauseidel.m如下：function [y,n]=gauseidel(A,b,x0,ep)if nargin==3ep=1.0e-6;elseif nargin<3errorreturnendD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1; %迭代次数while norm(y-x0)>=epx0=y;y=G*x0+f;n=n+1;end例8-6 用Gauss-Serdel迭代法求解例8-5中的线性方程组。
第8章 常微分方程的数值求解所谓的常微分方程就是把自变量t 和它的函数y 以及它的微商dy/dt 、d 2y/dt 2、…d n y/dt n 相联系的一个关系式0,...,,,22=⎪⎪⎭⎫ ⎝⎛n n dt y d dt y d dt dy y t f (1) 一个微分方程不只有一个或几个解，而是有无数个一簇解。
如一阶常微分方程dy/dt=1-e -t 的解为y(t)=t+e -t +C 。
C 为积分常数，C 取任意数值时，函数y(t)都满足微分方程。
因此，解有无数个。
如图在实际应用中，并不要求把所有的解都求出来，而是求满足某种指定条件的解。
这个条件通常称定解条件。
一个最重要的定解条件是初值条件。
对于上述方程，初值条件是：()00y t y =，()'00y dt t dy =，())2(0202y dt t dy =，…，())1(0101---=n n n y dt t dy (2) x 0是自变量的某个指定的“初值”，而0y 、'0y 、)2(0y 、…、)1(0-n y 则是未知函数及其到n －1阶微商的指定“初值”。
求解满足这样初值条件的微分方程问题称为初值问题。
如果能从方程(1)将n n dt y d 解出，则微分方程变成下面的形式⎪⎪⎭⎫ ⎝⎛=--1122,...,,,n n n n dt y d dt y d dt dy y t f dt y d (3) 这里的f 与式(1)中的f 不同，它是n+1个自变量的已知函数。
这种微分方程称为正规形微分方程。
而式(1)有时称为隐微分方程。
我们只考虑正规形微分方程，而且是一阶常微分方程的初值问题： ()y t f dtdy,=，()00y t y =， (4) 设函数()y t f ,在区域T t t ≤≤0，∞<u 内连续，并且存在常数L(Lipschitz 常数)，对所有],[00T t t ∈和y 1、y 2，有()()2121,,y y L y t f y t f -≤-，由常微分方程理论得知，初值问题(4)在区间[t 0, T]有唯一解，且连续可微。
微分方程的数值解法是一种离散化方法。
在这种方法中，事先取定自变量t 的离散值t 1，t 2，…，t N ，t j 称为节点，通常取成等距的，即t 1=t 0+h ，t 2=t 0+2h ，…，t N =t 0+Nh ，其中h>0称为步长，必要时，可以改变它的大小。
在这些节点上求出未知函数y(t)对应的值y(t 1)、y(t 2)、…y(t N )的近似值y 1、y 2、…y N 。
这些值称为初值问题的一个数值解。
常微分方程的数值求解方法很多，常用的有欧拉(Euler)方法、休恩(Heun)方法、泰勒(Taylor)级数法、龙格－库塔(Runge-Kutta)法、预测－校正(Predictor- Corrector)法。
欧拉(Euler)方法根据泰勒(Taylor)定理，如果函数y(t)、y ’(t)、y ’’(t)是连续的，在t=t 0附近可以展开为：()()()()()()2'''201000t t c y t t t y t y t y -+-+= ],[01t t c ∈ 令01t t h -=，并将()()()000,'t y t f t y =代入上式，得到()()()()()2'',210001h c y t y t hf t y t y ++=若步长取得足够小，可以忽略包含h 2的二次项()()()0001,t y t hf t y y +=——欧拉逼近将区间[a, b]分成M 个等距的子区间，子区间的长度即为步长h=(b-a)/M ，t k =a+kh, k=1, 2, …., M 或 t k+1= t k +h 。
重复上述欧拉逼近过程：()1112,y t hf y y +=；()2223,y t hf y y +=；h t t k k +=+1，()k k k k y t hf y y ,1+=+这些点就是对曲线y=y(t)的逼近。
在几何上，()()dt dy y t f =,表示曲线)(t y y =的斜率。
对于方程解的所有曲线，()y t f ,构成一个斜率场。
一个解曲线应该受斜率约束的。
沿着一条曲线运动，要从初始点(t 0, y 0)开始，考察在斜率场中该点的方向，沿着该方向移动一段距离，相当于在t 方向移动h 距离，在y 方向移动()00,y t hf ，这时移动到下一点(t 1, y 1)。
在(t 1, y 1)点，重复上述步骤，得到(t 2, y 2)，….，最后得到(t N , y N )。
当然，除了(t 0, y 0)外，其它点并不在精确的解曲线上，因此(t 1, y 1)，….，(t N , y N )是精确解曲线的近似值。
如果步长取得越短，近似值就逼近精确值。
欧拉方法的末端累积误差为()()()()12(2h O h c y a b =- 以下是在[0, 3]内，方程为()2'y t y -=、初值条件为()10=y 的不同步长欧拉解的比较。
从图中可以看出，步长取得越小，越逼近精确解。
休恩(Heun)方法对微分方程()y t f dt dy ,=从t 0到t 1积分，利用初始条件()00y t y =，得到()()()()()01110',t y t y dt t y dt t y t f t t t t -==⎰⎰→()()()()⎰+=1,01t t dt t y t f t y t y上式中的积分可以利用数值积分的梯形公式(()()[]b f a f ab S +-=21)表达。
现在取01t t h -=，应用梯形公式，得到()()()()()()()110001,,2t y t f t y t f ht y t y ++≈ (Heun-1) 注意，这时公式的右端也包含待定值y(t 1)。
可以用一个估计值来代替右边的y(t 1)。
这里采用欧拉方法得到这个估计值，()()()()0001,t y t hf t y t y +≈。
然后代入(Heun-1)式中，得到(t 1, y 1)。
()()()()()()00010001,,,2y t hf y t f t y t f ht y y +++=重复上述步骤，就得到一系列的点逼近解函数的精确值。
这种方法就是休恩(Heun)方法。
在这个方法中，用欧拉方法作为预报，然后用梯形公式进行校正。
它的一般步骤是()k k k k y t hf y p ,1+=+，h t t k k +=+1()()()1111,,2+++++=k k k k k k p t f y t f hy y 在几何上，用休恩方法得到的近似值比用欧拉方法得到的更接近准确值。
休恩方法的累积误差为()()()()222(12h O h c y a b =-以下是在[0, 2]内，方程为()2'y t y -=、初值条件为()10=y 的不同步长休恩方法解的比较。
可以看出，用休恩方法要比用欧拉方法收敛快。
泰勒(Taylor)级数法根据泰勒定理，如果函数y(t)及其直到N+1阶导数在区间[t 0, b]连续(即()],[01b t C t y N +∈)，则y(t)在],[0b t t t k ∈=处展开的N 次泰勒级数为()()()()()1,+++=+N k k N k k h O t y t hT t y h t y 其中，()()()()∑=-=Nj j k j k k N h j t y t y t T 11!,()()()()()t y t f t y j j ,1-=表示函数f 对t 的(j-1)次全导数：()f t y ='，()f f f t y y t +=''，()()()f f f f f f f f f t y y t y yy ty tt ++++=232，….()()()()()t y t f P t y N N ,1-=，⎪⎪⎭⎫⎝⎛∂∂+∂∂=y f t P 在每个子区间[t k , t k+1]应用泰勒定理就可以导出微分方程()y t f y ,'=的近似数值解。
这就是泰勒级数法。
N 次泰勒级数法的一般步骤是：!!3!2332211N h d h d h d h d y y N N k k +++++=+在k=0, 1, …, M-1的每一步，()()k j j t y d =，j=1, 2, …, N很显然，N 次泰勒级数法的误差量级是()1+N h O 。
因此，可以选择所需要的N ，使误差足够小。
如果N 固定，则可以找出满足要求的步长h 。
龙格－库塔(Runge-Kutta)法泰勒级数法的显著优点是误差是()1+N h O 量级的，而且N 可以选得较大使得误差较小。
但是，泰勒级数法的缺点是，要事先确定级数次数N 和计算高阶导数，通常这是很复杂的。
利用龙格－库塔法，即可以达到较高的精度，又可以避免高阶导数的计算。
由常微分方程理论知，初值问题(()y t f dt dy ,=，()00y t y =)可以写成等价的积分形式()()()⎰+=tt d y f y t y 0,0τττ。
于是，()()()()⎰++=+h t td y f t y h t y τττ,。
由积分中值定理知，()()()()h t y h t hf dt t y t f n n t t n nθθ++=⎰+,,1，10<<θ。
可得()()()()h t y h t hf t y h t y n n n n θθ+++=+,但是，()()h t y h t f n n θθ++,的值是无法计算的。
现在在区间],[h t t n n +取若干个点(例如s 个点)，用f 在这些点的值的线性组合来近似()()h t y h t f n n θθ++,的值，并具有尽可能高的精度。
具体地讲，就是用下面的公式代替上式()∑=++=si i i n n K W h t y y 11()()()⎪⎪⎩⎪⎪⎨⎧=⎪⎪⎭⎫ ⎝⎛++==∑-=s i K h t y h t f K t y t f K i j j ij n i ni n n ,....,2,,11βα 其中，W i 、αi 、βij 都是待定系数。
上两式称为s 级龙格－库塔法。
W i 、αi 、βij 的选择应该使误差具有尽可能高的阶。
通过龙格－库塔公式与泰勒级数公式比较，可以确定出待定系数W i 、αi 、βij 。
以s=4为例，这时关于W i 、αi 、βij 方程共有11个， 未知数有13个。
