MATLAB求解常微分方程数值解
利用MATLAB求解常微分方程数值解
目录
1. 内容简介 (1)
2. Euler Method(欧拉法)求解 (1)
. 显式Euler法和隐式Euler法 (2)
. 梯形公式和改进Euler法 (3)
. Euler法实用性 (5)
3. Runge-Kutta Method(龙格库塔法)求解 (6)
. Runge-Kutta基本原理 (6)
. MATLAB中使用Runge-Kutta法的函数 (8)
4. 使用MATLAB求解常微分方程 (8)
. 使用ode45函数求解非刚性常微分方程 (8)
. 刚性常微分方程 (9)
5. 总结 (10)
参考文献 (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程序见附录，能够得出精确解和数值解的图像，如图所示。
图显式Euler法精确解和数值解图像
从图中可以看出，显式Euler法在斜率很大的时候存在非常大的误差。本质上是Euler 法只计算了每一步差值中的一阶部分，由Taylor级数可知：
当公式中的二阶导数较大时就会产生明显的偏差，同时迭代过程中由于使用到上一部的结果，误差会在迭代中传播，因此这种Euler法在实际中是无法使用的，但是却给求解微分方程数值解提供了好的开始。
另外一种Euler法是隐式Euler法，其迭代公式是，它并没有解决上面所说的问题，同时它的计算更加繁琐，对于无法化简成显示迭代的公式时还需要用迭代法求解非线性方程。
为了解决上面的方法，就需要提高迭代公式中计算差值的阶数，下面介绍了梯形法和改进Euler法，它们都是二阶方法。
2.2.梯形公式和改进Euler法
梯形公式以及改进Euler法都属于二阶方法，下面证明它是二阶方法，使用两次Taylor 公式，将和展开：
将得到
从上式可以看出，梯形法的局部截断误差的主要部分是，是关于步长的三次式，这说明了梯形法取到了差值中的二次项，因此梯形法是二阶方法。
从上面可以得到梯形的迭代公式：
但是上式并不容易计算，因为上式中的为带求量，当无法化成显式形
式时，需要对上式进行迭代求解。因此梯形公式不易通过计算机编程求解，实际上改进的Euler法更容易求解。
改进Euler法迭代公式先通过显式Euler法求出一个估计值，通过这个估计值来计
算，然后方程就变成了显式方程，从而可以得到修正值，改进Euler法
更适合计算机编写程序，同样解决初值问题，详细MATLAB程序见附录2，得到的对比
图像如图所示。
图改进Euler法精确解与数值解对比
由于改进Euler法用来求解的并不是精确解，所以得到的导数会有一定误
差，因此改进Euler法的实际局部截断误差不仅仅是。
2.3.Euler法实用性
从图可以看出来一阶方法精确度非常差，基本上是无法用到实际工程中的，因此显式和隐式Euler法只是提供一种对微分方程求解的思想。
从图中得到的数值解相对图已经有了明显的改善，但是对于精度要求较高的工程问题，梯形法和改进Euler法这样的二阶方法同样是不满足要求的。但通过这两个方法，了解到提高方法取差值中的更高阶数项能够达到提高精度的目的，后面内容中的Runge-Kutta法就是由此思想而来。
将细节部分放大更能够看出两种方法的精度，如图所示(左图是Euler法，右图是改进Euler法)。
图两种方法细节部分
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)，得到的图像如图所示，图中解析解和数值解基本上重叠了，并且在处的误差已经小于了，而显式Euler法的误差为数量级，二阶改进Euler法误差在数量级。
图四阶四级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)如图所示，其精度与上面自
己编写的程序差不多，在处的误差小于。
图使用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 = ;
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)+,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 = ;
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)+,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 = ;
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.使用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::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 ]);