MATLAB求解常微分方程数值解

利用MATLAB求解常微分方程数值解

目录

1. 内容简介 (1)

2. Euler Method(欧拉法)求解 (1)

2.1. 显式Euler法和隐式Euler法 (2)

2.2. 梯形公式和改进Euler法 (3)

2.3. Euler法实用性 (4)

3. Runge-Kutta Method(龙格库塔法)求解 (5)

3.1. Runge-Kutta基本原理 (5)

3.2. MATLAB中使用Runge-Kutta法的函数 (7)

4. 使用MATLAB求解常微分方程 (7)

4.1. 使用ode45函数求解非刚性常微分方程 (8)

4.2. 刚性常微分方程 (9)

5. 总结 (9)

参考文献 (11)

附录 (12)

1. 显式Euler法数值求解 (12)

2. 改进Euler法数值求解 (12)

3. 四阶四级Runge-Kutta法数值求解 (13)

4.使用ode45求解 (14)

1.内容简介

把《高等工程数学》看了一遍，增加对数学内容的了解，对其中数值解法比较感兴趣，这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。理解模型然后列出微分方程，却对着方程无从下手，无法得出精确结果实在是让人难受的一件事情。

实际问题中更多遇到的是利用数值法求解偏微分方程问题，但考虑到先从常微分方程下手更为简单有效率，所以本文只研究常微分方程的数值解法。把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思，因此文章中不追求非常严格地证明，而是偏向如何利用工具实际求解出常微分方程的数值解，力求将课程上所学的知识真正地运用到实际方程的求解中去，在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。

文中求解过程中用到MATLAB进行数值求解，主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的，对各种方法进行MATLAB编程求解，并将求得的数值解与精确解对比，其中源程序在附录中。最后考察MATLAB中各个函数的适用范围，当遇到实际工程问题时能够正确地得到问题的数值解。

2.Euler Method(欧拉法)求解

Euler法求解常微分方程主要包括3种形式，即显式Euler法、隐式Euler法、梯形公式法，本节内容分别介绍这3种方法的具体内容，并在最后对3种方法精度进行对比，讨论Euler法的实用性。

本节考虑实际初值问题

使用解析法，对方程两边同乘以得到下式

两边同时求积分并采用分部积分得到解析解：

本节后面将对此方程进行求解，并与精确解进行对比，分析Euler的可行性。

2.1.显式Euler法和隐式Euler法

显式和隐式Euler法都属于一阶方法，显式Euler法的迭代公式简单，如下所示：

对过上述公式对式进行迭代，其中步长，计算之间的数值，迭代求解

的MATLAB程序见附录，能够得出精确解和数值解的图像，如图所示。

图2.1 显式Euler法精确解和数值解图像

从图2.1中可以看出，显式Euler法在斜率很大的时候存在非常大的误差。本质上是Euler 法只计算了每一步差值中的一阶部分，由Taylor级数可知：

当公式中的二阶导数较大时就会产生明显的偏差，同时迭代过程中由于使用到上一部的结果，误差会在迭代中传播，因此这种Euler法在实际中是无法使用的，但是却给求解微分方程数值解提供了好的开始。

另外一种Euler法是隐式Euler法，其迭代公式是，它并没有解决上面所说的问题，同时它的计算更加繁琐，对于无法化简成显示迭代的公式时还需要用迭代法求解非线性方程。

为了解决上面的方法，就需要提高迭代公式中计算差值的阶数，下面介绍了梯形法和改进Euler法，它们都是二阶方法。

2.2.梯形公式和改进Euler法

梯形公式以及改进Euler法都属于二阶方法，下面证明它是二阶方法，使用两次Taylor 公式，将和展开：

将得到

从上式可以看出，梯形法的局部截断误差的主要部分是，是关于步长的三次式，这说明了梯形法取到了差值中的二次项，因此梯形法是二阶方法。

从上面可以得到梯形的迭代公式：

但是上式并不容易计算，因为上式中的为带求量，当无法化成显式形

式时，需要对上式进行迭代求解。因此梯形公式不易通过计算机编程求解，实际上改进的Euler法更容易求解。

改进Euler法迭代公式先通过显式Euler法求出一个估计值，通过这个估计值来计算

，然后方程就变成了显式方程，从而可以得到修正值，改进Euler法更

适合计算机编写程序，同样解决初值问题，详细MATLAB程序见附录2，得到的对比图

像如图2.2所示。

图2.2 改进Euler法精确解与数值解对比

由于改进Euler法用来求解的并不是精确解，所以得到的导数会有一定误

差，因此改进Euler法的实际局部截断误差不仅仅是。

2.3.Euler法实用性

从图2.1可以看出来一阶方法精确度非常差，基本上是无法用到实际工程中的，因此显式和隐式Euler法只是提供一种对微分方程求解的思想。

从图2.2中得到的数值解相对图2.1已经有了明显的改善，但是对于精度要求较高的工程问题，梯形法和改进Euler法这样的二阶方法同样是不满足要求的。但通过这两个方法，了解到提高方法取差值中的更高阶数项能够达到提高精度的目的，后面内容中的Runge-Kutta法就是由此思想而来。

将细节部分放大更能够看出两种方法的精度，如图2.3所示(左图是Euler法，右图是改进Euler法)。

图2.3 两种方法细节部分

3.Runge-Kutta Method(龙格库塔法)求解

这一节的目的是弄清楚Runge-Kutta法(后面简称R-K法)的基本原理，并弄清楚MATLAB 中ode系列适合解决哪一种工程问题，希望达到的目的对于任何一个工程上提出的常微分方程问题都能够利用MATLAB求得误差在能够接受范围内的数值解。

3.1.Runge-Kutta基本原理

MATLAB中使用R-K法的函数有ode23和ode45，其中ode是ordinary differential equation 的缩写，其中ode23表示使用的是2阶3级R-K法，ode45使用的是4阶5级R-K法。

将局部截断误差表示如下式：

将上式中的方法称作阶级Runge-Kutta法。其中：

下面将构造出四阶四级经典R-K法，并通过MATLAB编程实现对初值问题求解。

四级四阶R-K法的本质是用四项级数来取得Taylor展开的前面四阶项，局部截断误差首项是步长的五次方，即：

从上述表达式中可以看出，实际上R-K法是希望通过线性组合来近似Taylor公式中的阶

导数组合，这样能够避免求解函数的阶导数，从而大量减小了计算量。实际计算中将用来表示，则利用二元Taylor公式展开，舍去，并将同次数项系数相等，得到经典四阶四级R-K法的公式：

上面从一阶到二阶的数值求法让误差有很大的改观，利用上述公式写出四阶四级的MATLABA程序(见附录3)，得到的图像如图3.1所示，图中解析解和数值解基本上重叠了，并且在处的误差已经小于了，而显式Euler法的误差为数量级，二阶改进Euler 法误差在数量级。

图3.1 四阶四级R-K法数值解和精确解对比

3.2.MATLAB中使用Runge-Kutta法的函数

MATLAB中利用R-K法求解常微分方程的函数主要是ode45，它用来求解非刚性常微分方程，是工程计算中首选函数。另外的ode23是2阶3级R-K法函数，它可以用来求解轻微刚性方程，其精度相对较低，但是效率较高。它们都是显式的R-K法，另外一个函数式ode23tb，它是隐式的R-K法可以用来求解非刚性的常微分方程。

4.使用MATLAB求解常微分方程

MATLAB中求解常微分方程共有7个函数，表4-1列举出它们的使用范围以及简介。

表4-1 MATLAB中求解常微分方程函数

通过表4-1了解到，当遇到一个常微分方程问题，首先使用ode45对其进行求解，如果求解非常慢或者无法求解，那么问题很可能是刚性问题，需要使用ode15s进行求解。对于精度要求非常高的问题可使用ode113尝试进行求解，对于其他的刚性问题，需要具体研究使用哪一个函数。

4.1.使用ode45函数求解非刚性常微分方程

ode45函数实现了变步长四阶五级Runge-Kutta-Felhberg算法(简称RKF法)。调用格式如下所示。

公式中是函数的句柄，是求解区间的起点和终点，是初始值，后面的

都是各种控制求解细节的选项。其中的可以是匿名函数、子函数、嵌套函数、单独M 文件等形式。实际上ode45能够求解常微分方程组，其中的使用状态变量表示方法，左侧是各个状态变量的导数，右侧是方程。ode45可以解决高阶常微分方程的数值求解问题，基本思想是将高阶通过变量替换变成低阶，然后列成下面所示的微分方程组，然后进行求解。

对于处置初值问题，使用MATLAB求解(代码见附录4)如图3.2所示，其精度与上

面自己编写的程序差不多，在处的误差小于。

图3.2 使用ode45求解细节

4.2.刚性常微分方程

如上所述，ode45函数能够求解的是非刚性问题，而刚性(stiff)问题使用ode45求解会产生求解缓慢，或者完全无法求解的情况，这种情况可以判断这个方程是刚性的。可以考虑使用ode23或者ode23s函数进行求解。

刚性常微分方程指的是其特征根相差很大，但似乎没有非常有效判断一个方程是刚性还是非刚性的办法，需要实际求解的时候通过不同的函数进行求解比较。

5.总结

文本主要研究了如何求常微分方程的数值解，从开始的一阶显式Euler法，到四阶R-K 法，精度提高了很多。同时通过编写MATLAB代码来对具体初值问题进行求解，对课本上将

的机制有了更加深入的了解，对具体算法过程有了深入的认识，从中学到了很多。

另外文中总结了MATLAB中求解常微分方程数值解函数的用法，在面对不同的工程问题的时候，怎样去挑选正确的函数来求解具体问题，更多的是学会了如何去解决实际问题。

实际问题中更多是由偏微分方程进行描述的，由于时间和个人数学水平原因，文中没有对偏微分方程数值求解进行研究，这是文章的不足，有待在后续的学习当中对其进行解决。

参考文献

[1]于寅.高等工程数学[M].武汉:华中科技大学出版社.2012

[2]李红.数值分析(第2版)[M].武汉:华中科技大学出版社.2010

[3]胡建伟,汤怀民.微分方程数值方法(第二版)[M].北京:科学出版社.2007

[4]薛定宇,陈阳泉.高等应用数学问题的MATLAB求解(第二版)[M].北京:清华大学出版社.2008

附录

1.显式Euler法数值求解

clear,clc %y1中存储0~5中的数值解y2中存储-5~0中的数值解dyx=@(x,y) (-2*y-4*x); %方程右边函数

fx=@(x) (exp(-2*x)-2*x+1); %解析解

y1(1)=2;y2(1)=2;

h = 0.1;

i=2;

for x=0:h:5-h

y1(i)=y1(i-1)+h*dyx(x,y1(i-1));

i =i +1;

end

i=2;

for x=0:-h:-5+h

y2(i) = y2(i-1)-h*dyx(x,y2(i-1));

i=i+1;

end

x=-5:h:5;

y=fx(x);

ye=[fliplr(y2(2:end)) 2 y1(2:end)];

plot(x,y,x,ye,'--');

legend('\fontsize{12}\fontname{楷体}精确解',...

'\fontsize{12}\fontname{楷体}数值解','interpreter','latex');

text(x(8)+0.1,y(8),'$\leftarrow e^{-2x}-2x+1$','interpreter','latex','fontsize',15);

2.改进Euler法数值求解

clear,clc %y1中存储0~5中的数值解y2中存储-5~0中的数值解dyx=@(x,y) (-2*y-4*x); %方程右边函数

fx=@(x) (exp(-2*x)-2*x+1); %解析解

y1(1)=2;y2(1)=2;

yy1(1)=2;yy2(1)=2; %yy1和yy2中存储的为估计值

h = 0.1;

i=2;

for x=0:h:5-h

yy1(i)=y1(i-1)+h*dyx(x,y1(i-1));

yy1(i)=dyx(x+h,yy1(i));

y1(i)=y1(i-1)+h/2*(dyx(x,y1(i-1))+yy1(i));

i =i +1;

end

i=2;

for x=0:-h:-5+h

yy2(i) = y2(i-1)-h*dyx(x,y2(i-1));

yy2(i) = dyx(x-h,yy2(i));

y2(i) =y2(i-1)-h/2*(dyx(x,y2(i-1))+yy2(i));

i=i+1;

end

x=-5:h:5;

y=fx(x);

ye=[fliplr(y2(2:end)) 2 y1(2:end)];

plot(x,y,x,ye,'--');

set(gcf,'color','white');

legend('\fontsize{12}\fontname{楷体}精确解',...

'\fontsize{12}\fontname{楷体}数值解','interpreter','latex');

text(x(8)+0.1,y(8),'$\leftarrow e^{-2x}-2x+1$','interpreter','latex','fontsize',15);

3.四阶四级Runge-Kutta法数值求解

clear,clc %y1中存储0~5中的数值解y2中存储-5~0中的数值解dyx=@(x,y) (-2*y-4*x); %方程右边函数,即f(x,y)

fx=@(x) (exp(-2*x)-2*x+1); %解析解

y1(1)=2;y2(1)=2;

h = 0.1;

i=2;

for x=0:h:5-h

k1=dyx(x,y1(i-1));

k2=dyx(x+h/2,y1(i-1)+h/2*k1);

k3=dyx(x+h/2,y1(i-1)+h/2*k2);

k4=dyx(x+h,y1(i-1)+h*k3);

y1(i)=y1(i-1)+h/6*(k1+2*k2+2*k3+k4);

i =i +1;

end

i=2;

for x=0:-h:-5+h

k1=dyx(x,y2(i-1));

k2=dyx(x-h/2,y2(i-1)-h/2*k1);

k3=dyx(x-h/2,y2(i-1)-h/2*k2);

k4=dyx(x-h,y2(i-1)-h*k3);

y2(i)=y2(i-1)-h/6*(k1+2*k2+2*k3+k4);

i=i+1;

end

x=-5:h:5;

y=fx(x);

ye=[fliplr(y2(2:end)) 2 y1(2:end)];

subplot(1,2,1);

plot(x,y,x,ye,'--');

set(gcf,'color','white');

legend('\fontsize{12}\fontname{楷体}精确解',...

'\fontsize{12}\fontname{楷体}数值解','interpreter','latex'); hold on;

subplot(1,2,2);

plot(x,y,x,ye,'--');

legend('\fontsize{12}\fontname{楷体}精确解',...

'\fontsize{12}\fontname{楷体}数值解','interpreter','latex'); axis([-5 -4.9995 2.2030e4 2.204e4]);

4.使用ode45求解

dyx=@(x,y) (-2*y-4*x); %方程右边函数

fx=@(x) (exp(-2*x)-2*x+1); %解析解

[xx,yy]=ode45(dyx,[0 5],2);

[xxx,yyy]=ode45(dyx,[0 -5],2),

x=-5:0.1:5;

y=fx(x);

sx=[flipud(xxx) xx];

sy=[flipud(yyy) yy];

subplot(1,2,1);

plot(x,y,sx,sy,'--');

subplot(1,2,2);

plot(x,y,sx,sy,'--');

set(gcf,'color','white');

axis([-5 -4.9997 2.203e4 2.205e4]);

用MATLAB解常微分方程

实验四 求微分方程的解 一、问题背景与实验目的 实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）． 对微分方程（组）的解析解法(精确解)，Matlab 有专门的函数可以用，本实验将作一定的介绍． 本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍 Euler 折线法． 二、相关函数（命令）及简介 1．dsolve ('equ1','equ2',…)：Matlab 求微分方程的解析解．equ1、equ2、…为方程（或条件）．写方程（或条件）时用 Dy 表示y 关于自变量的一阶导数，用用 D2y 表示 y 关于自变量的二阶导数，依此类推． 2．simplify(s )：对表达式 s 使用 maple 的化简规则进行化简． 例如： syms x simplify(sin(x)^2 + cos(x)^2) ans=1 3．[r,how]=simple(s)：由于 Matlab 提供了多种化简规则，simple 命令就是对表达式 s 用各种规则进行化简，然后用 r 返回最简形式，how 返回形成这种形式所用的规则． 例如： syms x [r,how]=simple(cos(x)^2-sin(x)^2) r = cos(2*x) how = combine 4．[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解． 说明： (1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一． (2) odefun 是显式常微分方程：?????==0 0)() ,(y t y y t f dt dy (3) 在积分区间 tspan =],[0f t t 上，从0t 到f t ，用初始条件0y 求解．

用Matlab解微分方程

用Matlab 软件求解微分方程 1．解析解 （1）一阶微分方程 求21y dx dy +=的通解：dsolve('Dy=1+y^2','x') 求y x dx dy -+=21的通解：dsolve('Dy=1+x^2-y','x') 求?????=+=1 )0(12y y dx dy 的特解：dsolve('Dy=1+y^2',’y(0)=1’,'x') （2）高阶微分方程 求解???-='==-+'+''. 2)2(,2)2(,0)(222πππy y y n x y x y x 其中，21=n ，命令为： dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x') 求042=-+'-'''x y y y 的通解，命令为： dsolve('D3y-2*Dy+y-4*x=0','x') 输出为： ans=8+4*x+C1*exp(x)+C2*exp(-1/2*(5^(1/2)+1)*x)+C3*exp(1/2*(5^(1/2)-1)*x) （3）一阶微分方程组 求???+-='+='). (3)(4)(),(4)(3)(x g x f x g x g x f x f 的通解：[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','x') 输出为： f =exp(3*x)*(cos(4*x)*C1+sin(4*x)*C2) g =-exp(3*x)*(sin(4*x)*C1-cos(4*x)*C2) 若再加上初始条件1)0(,0)0(==g f ，则求特解： [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1','x') 输出为： f =exp(3*x)*sin(4*x) g =exp(3*x)*cos(4*x) 2．数值解 （1）一阶微分方程

Matlab求解微分方程(组)及偏微分方程(组)

第四讲Matlab求解微分方程(组) 理论介绍:Matlab求解微分方程(组)命令 求解实例:Matlab求解微分方程(组)实例 实际应用问题通过数学建模所归纳得到得方程,绝大多数都就是微分方程,真正能得到代数方程得机会很少、另一方面,能够求解得微分方程也就是十分有限得,特别就是高阶方程与偏微分方程(组)、这就要求我们必须研究微分方程(组)得解法:解析解法与数值解法、 一.相关函数、命令及简介 1、在Matlab中,用大写字母D表示导数,Dy表示y关于自变量得一阶导数,D2y 表示y关于自变量得二阶导数,依此类推、函数dsolve用来解决常微分方程(组)得求解问题,调用格式为: X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解、 注意,系统缺省得自变量为t 2、函数dsolve求解得就是常微分方程得精确解法,也称为常微分方程得符号解、但就是,有大量得常微分方程虽然从理论上讲,其解就是存在得,但我们却无法求出其解析解,此时,我们需要寻求方程得数值解,在求常微分方程数值解方 面,MATLAB具有丰富得函数,我们将其统称为solver,其一般格式为: [T,Y]=solver(odefun,tspan,y0) 说明:(1)solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb、ode15i之一、 (2)odefun就是显示微分方程在积分区间tspan上从到用初始条件求解、 (3)如果要获得微分方程问题在其她指定时间点上得解,则令tspan(要求就是单调得)、 (4)因为没有一种算法可以有效得解决所有得ODE问题,为此,Matlab提供了多种求解器solver,对于不同得ODE问题,采用不同得solver、 表1 Matlab中文本文件读写函数

用Matlab软件求常微分方程的解或通解

《高等数学》实验报告 实验人员:系(班): 学号: 姓名: 实验地点:电教楼五号机房 实验名称:Matlab 高等数学实验 实验时间:2014-6-3 16:30--18:30 实验名称:用Matlab 软件求常微分方程的解(或通解) 实验目的:熟练掌握Matlab 软件求常微分方程的解(或通解) 实验内容:(给出实验程序与运行结果) 1、求微分方程的特解. 1、?? ?? ?===+-10)0(,6)0(034'2 2y y y dx dy dx y d 程序:>> dsolve('D2y-4*Dy+3*y','y(0)=6,Dy(0)=10','x') ans = 4*exp(x)+2*exp(3*x) 吕梁学院《高等数学》实验报告 情况试高中

2、?? ???===++0)0(,2)0(044'2 2y y y dx dy dx y d 程序:>>dsolve('4*D2y+4*Dy+y','y(0)=2,Dy(0)=0','x') ans = 2*exp(-1/2*x)+exp(-1/2*x)*x 3、?? ???===++15)0(',0)0(029422y y y dx dy dx y d 程序:>>dsolve('D2y+4*Dy+29*y=0','y(0)=9,Dy(0)=15','x') ans = 33/5*exp(-2*x)*sin(5*x)+9*exp(-2*x)*cos(5*x) 4、?? ???===+-3)0(',0)0(013422y y y dx dy dx y d 程序:>>dsolve('D2y-4*dy+13*y=0','y(0)=0','Dy(0)=3','x') ans = 3/13*sin(13^(1/2)*x)*13^(1/2)-4/13*cos(13^(1/2)*x)*dy+4/13*dy 5、?? ???-===--5)0(',0)0(04322y y y dx dy dx y d 程序:>>dsolve('D2y-3*Dy-4*y','y(0)=0,Dy(0)=-5','x') ans = exp(-x)-exp(4*x)

MATLAB求解常微分方程数值解

利用MATLAB求解常微分方程数值解

目录 1. 内容简介 (1) 2. Euler Method(欧拉法)求解 (1) 2.1. 显式Euler法和隐式Euler法 (2) 2.2. 梯形公式和改进Euler法 (3) 2.3. Euler法实用性 (4) 3. Runge-Kutta Method(龙格库塔法)求解 (5) 3.1. Runge-Kutta基本原理 (5) 3.2. MATLAB中使用Runge-Kutta法的函数 (7) 4. 使用MATLAB求解常微分方程 (7) 4.1. 使用ode45函数求解非刚性常微分方程 (8) 4.2. 刚性常微分方程 (9) 5. 总结 (9) 参考文献 (11) 附录 (12) 1. 显式Euler法数值求解 (12) 2. 改进Euler法数值求解 (12) 3. 四阶四级Runge-Kutta法数值求解 (13) 4.使用ode45求解 (14)

1.内容简介 把《高等工程数学》看了一遍，增加对数学内容的了解，对其中数值解法比较感兴趣，这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。理解模型然后列出微分方程，却对着方程无从下手，无法得出精确结果实在是让人难受的一件事情。 实际问题中更多遇到的是利用数值法求解偏微分方程问题，但考虑到先从常微分方程下手更为简单有效率，所以本文只研究常微分方程的数值解法。把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思，因此文章中不追求非常严格地证明，而是偏向如何利用工具实际求解出常微分方程的数值解，力求将课程上所学的知识真正地运用到实际方程的求解中去，在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。 文中求解过程中用到MATLAB进行数值求解，主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的，对各种方法进行MATLAB编程求解，并将求得的数值解与精确解对比，其中源程序在附录中。最后考察MATLAB中各个函数的适用范围，当遇到实际工程问题时能够正确地得到问题的数值解。 2.Euler Method(欧拉法)求解 Euler法求解常微分方程主要包括3种形式，即显式Euler法、隐式Euler法、梯形公式法，本节内容分别介绍这3种方法的具体内容，并在最后对3种方法精度进行对比，讨论Euler法的实用性。 本节考虑实际初值问题 使用解析法，对方程两边同乘以得到下式

Matlab求解微分方程(组)及偏微分方程(组)

第四讲 Matlab 求解微分方程（组） 理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介 1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解. 注意，系统缺省的自变量为t 2.函数dsolve 求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为： [T,Y]=solver(odefun,tspan,y0) 说明：(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一. (2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,, ,f t t t t 上的解，则令 tspan 012[,,,]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供

(完整版)实验七用matlab求解常微分方程

实验七 用matlab 求解常微分方程 一、实验目的： 1、熟悉常微分方程的求解方法，了解状态方程的概念； 2、能熟练使用dsolve 函数求常微分方程（组）的解析解； 3、能熟练应用ode45\ode15s 函数分别求常微分方程的非刚性、刚性的数值解； 4、掌握绘制相图的方法 二、预备知识： 1．微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为 0),,",',,()(=n y y y y t F Λ 如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为 ) ()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++--Λ 若上式中的系数 n i t a i ,,2,1),(Λ=均与t 无关，称之为常系数。 2．常微分方程的解析解 有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dt dy 可化为 dt y dy =+1,两边积分可得通解为 1-=t ce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。 一阶常微分方程与高阶微分方程可以互化，已给一个n 阶方程 ),,",',()1()(-=n n y y y t f y Λ 设 ) 1(21,,',-===n n y y y y y y Λ，可将上式化为一阶方程组 ?????????====-),,,,(''''2113221n n n n y y y t f y y y y y y y ΛΛ 反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。 3．微分方程的数值解法 除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题

用Matlab件求常微分方程解(或通解)

用Matlab件求常微分方程解(或通解)

————————————————————————————————作者：————————————————————————————————日期：

《高等数学》实验报告 实验人员:系(班): 学号: 姓名: 实验地点:电教楼五号机房 实验名称:Matlab 高等数学实验 实验时间:2014-6-3 16:30--18:30 实验名称:用Matlab 软件求常微分方程的解(或通解) 实验目的:熟练掌握Matlab 软件求常微分方程的解(或通解) 实验内容:(给出实验程序与运行结果) 一、求微分方程的特解. 1、?? ???===+-10)0(,6)0(034'22y y y dx dy dx y d 程序:>> dsolve('D2y-4*Dy+3*y','y(0)=6,Dy(0)=10','x') ans = 4*exp(x)+2*exp(3*x) 吕梁学院《高等数学》实验报告

2、?? ???===++0)0(,2)0(044'22y y y dx dy dx y d 程序:>>dsolve('4*D2y+4*Dy+y','y(0)=2,Dy(0)=0','x') ans = 2*exp(-1/2*x)+exp(-1/2*x)*x 3、?? ???===++15)0(',0)0(029422y y y dx dy dx y d 程序:>>dsolve('D2y+4*Dy+29*y=0','y(0)=9,Dy(0)=15','x') ans = 33/5*exp(-2*x)*sin(5*x)+9*exp(-2*x)*cos(5*x) 4、?? ???===+-3)0(',0)0(013422y y y dx dy dx y d 程序:>>dsolve('D2y-4*dy+13*y=0','y(0)=0','Dy(0)=3','x') ans = 3/13*sin(13^(1/2)*x)*13^(1/2)-4/13*cos(13^(1/2)*x)*dy+4/13*dy 5、?? ???-===--5)0(',0)0(04322y y y dx dy dx y d 程序:>>dsolve('D2y-3*Dy-4*y','y(0)=0,Dy(0)=-5','x') ans = exp(-x)-exp(4*x)

matlab求解常微分方程.docx

用matlab求解常微分方程 在MATLAB中，由函数dsolve()解决常微分方程(组)的求解问题，其具体格式如下：r 二dsolve('eql,eq2,???字condl,cond2,?.?；V) 匕ql,eq2,???*为微分方程或微分方程组，，condl,cond2,.??；是初始条件或边界条件，P是独立变量，默认的独立变量是讥 函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。 dy _1 例1：求解常微分方程莎一石的MATLAB程序为： dsolve(* Dy=l/(x+y) 1r!x1), 注意，系统缺省的自变量为t,因此这里要把自变量写明。其中：Y=lambertw(X)表示函数关系Y*exp(Y)二X。 例2：求解常微分方程E'-y— 0的MATLAB程序为： Y2=dsolve(1y*D2y-Dy A2=01, 1x f) Y2=dsolve(!D2y*y-Dy A2=0 J )

我们看到有两个解，其中一个是常数0。 dx 心 ？ —+ 5x + y = e dt 空_兀_3『= g2f 例3：求常微分方程组〔力 ' 通解的MATLAB 程序为： [X,Y]=dsolve(f Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t) 1, 111 ) [X,Y]=dsolve(1 Dx+2 *x-Dy=l0 * cos(t),Dx+Dy+2 *y=4 *exp(- 2*t) T ,f x(0)=2,y(0)=0f ,f t T ) 以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。但是，我们知 道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析 解，此吋，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰 富的函数，我们将其统称为solver,其一般格式为： i°cosr, 7 =2 例4：求常微分方程组 y = 0 z 通解的MATLAB 程序为:

MATLABEuler法解常微分方程

Euler法解常微分方程 Euler法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2计算判断是否成立，成立转到Step 3，否则继续进行Step 4 Step 3 计算 Step 4 Euler法解常微分方程算程序： function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x

Step 3 （1）做显性Euler预测 （2）将带入 Step 4计算判断是否成立，成立返回Step 3，否则继续进行Step 5 Step 5 改进Euler法解常微分方程算程序： function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x' y=y' 例：用改进Euler法计算下列初值问题（取步长h=0.25） 输入：fun=inline('-x*y^2') gaijineuler2(fun,2,[0 5],0.25) 得到： x = 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500

数学应用软件作业6-用Matlab求解微分方程(组)的解析解和数值解

数学应用软件作业6-用Matlab 求解微分方程(组)的解析解和数值解

注：上机作业文件夹以自己的班级姓名学号命名，文件夹包括如下上机报告和Matlab程序。 上机报告模板如下： 佛山科学技术学院 上机报告 课程名称数学应用软件 上机项目用Matlab求解微分方程(组)的解析解和数值解 专业班级姓名学号 一. 上机目的 1.了解求微分方程（组）的解的知识。 2.学习Matlab中求微分方程的各种解的函数，如 dsolve命令、ode45函数等等，其中注意把方程化为新的方程的形式。 3.掌握用matlab编写程序解决求解微分方程的问 题。 二. 上机内容 1、求高阶线性齐次方程：y’’’-y’’-3y’+2y=0。 2、求常微分方程组

2 210cos,2 24,0 t t t dx dy x t x dt dt dx dy y e y dt dt = - = ? +-== ?? ? ?++== ?? 3、求解 分别用函数ode45和ode15s计算求解，分别画出图形，图形分别标注标题。 4、求解微分方程 ,1 )0( ,1 '= + + - =y t y y 先求解析解，在[0,1]上作图； 再用ode45求数值解(作图的图形用“o”表示)，在同一副图中作图进行比较，用不同的颜色表示。 三. 上机方法与步骤 给出相应的问题分析及求解方法，并写出Matlab 程序，并有上机程序显示截图。 题1：直接用命令dsolve求解出微分方程的通解。 Matlab程序：

dsolve('D3y-D2y-3*Dy+2*y','x') 题2：将微分方程组改写为 5cos2exp(2) 5cos2exp(2) (0)2,(0)0 dx t t x y xt dy t t x y dt x y ? =+--- ? ? ? =-+-+- ? ? == ? ? ? ， 再用命令dsolve求解微分方程的通解。 Matlab程序： 建立timu2.m如下： [x,y]=dsolve('Dx=5*cos(t)+2*exp(-2*t)-x-y','Dy=-5*cos(t)+2*exp(-2*t)+x-y ','x(0)=2,y(0)=0','t') x=simple(x) y=simple(y)

matlab简介(解常微分方程绘制函数图像)

MATLAB简介 MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。 一、基本功能 MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。 二、特点 1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来； 2) 具有完备的图形处理功能,实现计算结果和编程的可视化; 3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握； 4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ,为用户提供了大量方便实用的处理工具。 三、优势 1.友好的工作平台编程环境 MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。 2.强大的科学计算机数据处理能力 MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能，可以用它来代替底层编程语言，如C和C++ 。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。

用matlab求解常微分方程

实验六 用matlab 求解常微分方程 1．微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为 0),,",',,()(=n y y y y t F 如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为 )()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++-- 若上式中的系数n i t a i ,,2,1),( =均与t 无关，称之为常系数。 2．常微分方程的解析解 有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dt dy 可化为 dt y dy =+1,两边积分可得通解为 1-=t ce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。 一阶常微分方程与高阶微分方程可以互化，已给一个n 阶方程 ),,",',()1()(-=n n y y y t f y 设)1(21,,',-===n n y y y y y y ，可将上式化为一阶方程组 ?????????====-),,,,(''''2113221n n n n y y y t f y y y y y y y 反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。 3．微分方程的数值解法 除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题 ???=<<=000)()),(,()('y t y t t t t y t f t y f

利用matlab编写S函数求解微分方程

利用matlab编写S函数求解微分方程自动化专业综合设计报告 自动化专业综合设计报告

函数求解微S编写设计题目：利用 matlab 分方程 自动化系统仿真实验室所在 实验室： 郭卫平 指导教师： 律迪迪学生姓名 200990519114 班级文自0921 学号 成绩评定： 自动化专业综合设计报告

一、设计目的 了解使用simulink的扩展工具——S-函数，s函数可以利用matlab的丰富资源，而不仅仅局限于simulink提供的模块，而用c或c++等语言写的s函数还可以实现对硬件端口的操作，还可以操作windows API 等的，它的魅力在于完美结合了simulink 框图简洁明快的特点和编程灵活方便的优点，提供了增强和扩展sinulink能力的强大机制，同时也是使用RTW实现实时仿真的关键。 二、设计要求 求解解微分方程 y'=y-2x/y 自动化专业综合设计报告 y(0)=1 要求利用matlab编写S函数求解 三、设计内容（可加附页） 【步骤1】获取状态空间表达式。

在matlab中输入 dsolve(‘Dy=y-2*x/y','y(0)=1'， 'x') 得到 y=(2*x+1).^(1/2); 【步骤2】建立s函数的m文件。 利用21·用S函数模板文件。 以下是修改之后的模板文件sfuntmpl.m 的内容。 function [sys,x0,str,ts] = sfuntmpl(t,x,u,flag) %SFUNTMPL S-function M-file General template define you can With % M-file S-functions, you own ordinary differential system equations (ODEs), discrete % equations, and/or just about any type of algorithm to be used within a %

数学模型之微分方程及其MATLAB求解

数学模型之微分方程及其MATLAB求解 ---卫星轨迹等经典例题求解分析1. 考虑初值问题画图 y'''?3y ''?y 'y = 0 y(0) = 0 y '(0) =1 y ' '(0) = ?1 2、 3、 【实验步骤与程序】 1. M -文件建立m函数文件

function y=f(t,x) y=[x(2);x(3);9*x(3)^2+x(1)*x(2)]; 求解微分方程，命令如下： x0=[0;1;-1]; [t,y]=ode45(@mm,[0,2.5],x0); plot(y(:,1),y(:,2)); figure(2); plot3(y(:,1),y(:,2),y(:,3))

2、M -文件建立m函数文件 function dx=appollo(t,x) mu=1/82.45; mustar=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mustar)^2+x(3)^2); dx=[x(2) 2*x(4)+x(1)-mustar*(x(1)+mu)/r1^3-mu*(x(1)-mustar)/r2^3 x(4) -2*x(2)+x(3)-mustar*x(3)/r1^3-mu*x(3)/r2^3];

求解微分方程，命令如下： x0=[1.2;0;0;-1.04935751]; options=odeset('reltol',1e-8); [t,y]=ode45(@appollo,[0,20],x0,options); plot(y(:,1),y(:,3)) title('Appollo卫星运动轨迹') xlabel('x') ylabel('y')

Matlab解微分方程(ODE+PDE)

常微分方程： 1 ODE解算器简介(ode**) 2 微分方程转换 3 刚性/非刚性问题(Stiff/Nonstiff) 4 隐式微分方程(IDE) 5 微分代数方程(DAE) 6 延迟微分方程(DDE) 7 边值问题(BVP) 偏微分方程(PDEs)Matlab解法 偏微分方程： 1 一般偏微分方程组(PDEs)的命令行求解 2 特殊偏微分方程(PDEs)的PDEtool求解 3 陆君安《偏微分方程的MATLAB解法 先来认识下常微分方程(ODE)初值问题解算器(solver) [T,Y,TE,YE,IE] = odesolver(odefun,tspan,y0,options) sxint = deval(sol,xint) Matlab中提供了以下解算器： 输入参数： odefun：微分方程的Matlab语言描述函数，必须是函数句柄或者字符串，必须写成Matlab

规范格式(也就是一阶显示微分方程组)，这个具体在后面讲解 tspan=[t0 tf]或者[t0,t1,…tf]：微分变量的范围，两者都是根据t0和tf的值自动选择步长，只是前者返回所有计算点的微分值，而后者只返回指定的点的微分值，一定要注意对于后者tspan必须严格单调，还有就是两者数据存储时使用的内存不同(明显前者多)，其它没有任何本质的区别 y0=[y(0),y’(0),y’’(0)…]：微分方程初值，依次输入所有状态变量的初值，什么是状态变量在后面有介绍 options：微分优化参数，是一个结构体，使用odeset可以设置具体参数，详细内容查看帮助 输出参数： T：时间列向量，也就是ode**计算微分方程的值的点 Y：二维数组，第i列表示第i个状态变量的值，行数与T一致 在求解ODE时，我们还会用到deval()函数，deval的作用就是通过结构体solution计算t 对应x值，和polyval之类的很相似！ 参数格式如下： sol：就是上次调用ode**函数得道的结构体解 xint：需要计算的点，可以是标量或者向量，但是必须在tspan范围内 该函数的好处就是如果我想知道t=t0时的y值，不需要重新使用ode计算，而直接使用上次计算的得道solution就可以 [教程] 微分方程转换为一阶显示微分方程组方法 好，上面我们把Matlab中的常微分方程(ODE)的解算器讲解的差不多了，下面我们就具体开始介绍如何使用上面的知识吧！ 现实总是残酷的，要得到就必须先付出，不可能所有的ODE一拿来就可以直接使用，因此，在使用ODE解算器之前，我们需要做的第一步，也是最重要的一步，借助状态变量将微分

matlabeuler法解常微分方程

Euler 法解常微分方程 Euler 法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2计算h n n +=判断b n ≤是否成立，成立转到Step 3，否则继续进行Step 4 Step 3 计算),(1n n n n y x hf y y +=+ Step 4 ),(1n n n n y x hf y y +=+ Euler 法解常微分方程算程序： function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x

指导教师： 年 月 日 改进Euelr 法解常微分方程 改进Euler 法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2 取一个以h 为步长，a ，b 分别为左右端点的矩阵 Step 3 （1）做显性Euler 预测),(1n n i i y x hf y y +=+ （2）将1+i y 带入)],(),([2 h 111+++++=i i i i i i y x f y x f y y Step 4计算h n n +=判断b n ≤是否成立，成立返回Step 3，否则继续进行Step 5 Step 5 )],(),([2 h 111+++++=i i i i i i y x f y x f y y 改进Euler 法解常微分方程算程序： function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x'

Matlab求解微分方程(组)及偏微分方程(组).

第四讲 Matlab 求解微分方程（组） 理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介 1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解. 注意，系统缺省的自变量为t 2.函数dsolve 求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为： [T,Y]=solver(odefun,tspan,y0) 说明：(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一. (2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,,,f t t t t 上的解，则令 tspan 012[,,, ]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供了多种求解器solver ，对于不同的ODE 问题，采用不同的solver.

Matlab求解常微分方程边值问题的方法

Matlab 求解常微分方程边值问题的方法：bvp4c 函数 常微分方程的边值问题，即boundary value problems ，简称BVP 问题，是指表达形式为 (,)((),())0'=??=?y f x y g y a y b 或(,,)((),(),)0'=??=? y f x y p g y a y b p 的方程组（p 是未知参数），在MA TLAB 中使用积分器bvp4c 来求解。 [命令函数] bvp4c [调用格式] sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) sol 为一结构体，sol.x 、sol.y 、sol.yp 分别是所选择的网格点及其对应的y(x)与y'(x)数值; bvp4c 为带边值条件常微分方程积分器的函数命令；odefun 为描述微分方程组的函数文件；bcfun 为计算边界条件g(f(a),f(b),p)=0的函数文件；solinit 为一结构体，solinit.x 与solinit.y 分别是初始网格的有序节点与初始估计值，边界值条件分别对应a=solinit.x(l)和b=solinit.x(end)； options 为bvpset 命令设定的可选函数，可采用系统默认值；p1, p2…为未知参数。 例 求常微分方程0''+=y y 在(0)2=y 与(4)2=-y 时的数值解。 [解题过程] 仍使用常用方法改变方程的形式： 令1=y y ，21'=y y ，则原方程等价于标准形式的方程组1221 ?'=??'=-??y y y y ； 将其写为函数文件twoode.m ； 同时写出边界条件函数对应文件twobc.m ； 分别使用结构solinit 和命令bvp4c 确定y-x 的关系； 作出y-x 的关系曲线图。 [算例代码] solinit =bvpinit(linspace(0,4,5),[1 0]); % linspace(0,4,5)为初始网格，[1,0]为初始估计值 sol=bvp4c(@twoode,@twobc,solinit); % twoode 与twobc 分别为微分方程与边界条件的函数，solinit 为结构 x=linspace(0,4); %确定x 范围 y=deval(sol,x); %确定y 范围 plot(x,y(1,:)); %画出y-x 的图形 %定义twoode 函数（下述代码另存为工作目录下的twoode.m 文件） function dydx= twoode(x,y) %微分方程函数的定义 dydx =[y(2) -abs(y(1))]; %定义twobc 函数（下述代码另存为工作目录下的twobc.m 文件） function res= twobc(ya,yb); %边界条件函数的定义 res=[ya(1);yb(1)+2];