080数学基础：常微分方程(组)数值解法.doc

数学基础： 微分方程（组）数值解法 ━基于ＭＡＴＬＡＢ的实现与分析 §1 常微分方程初值问题的数值解法 微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton第二定律： 　　　　　　　　　　（１） 和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题： 　　　　　　　　　　　　　（２） 的解为： 　　　　　　　　　　　（３） 但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实： １、绝大多数的常微分方程和常微分方程组得不到解析解； ２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）； 如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。 一般的一阶常微分方程初值问题是指如下的一阶常微分方程的定解问题： 　　　　　　　　　（７） 微分方程（组）的初值问题通常是对一动态过程演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。 一、最简单的数值解法——Euler 方法 假设要求在点（时刻），，处初值问题（７）的解的近似值。首先对式（７）的两端积分，得 　　　　　　　　　（８） 对于（８）的右边，如果被积函数用积分下限处的函数值代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面 积），则有 　　　　　　　　（９） 进而得到下式给出的递推算法—Euler 方法 　　　　　　　　（１０） 例１　用Euler 方法解如下初值问题，取， 　　　　　　　　　　　　 解：由（１０）得 　　　　　　　　　　　　　 结果如下： Euler_Method 如果取，其结果如下图所示： Euler_Method 附　应用程序　Euler_Method.m %Euler Method for the Example: y = y - 2sin(t)y^3; y(0)=1; clear all [T,z]=ode45(Euler_Method_file,[0;3],1); t=0:0.3:3; %t=0:0.1:3; %t=0:0.3:3; Lt=length(t); y=zeros(1,Lt); y(1)=1; for k=2:Lt y(k)=1.3*y(k-1)-0.6*sin(t(k-1))*y(k-1)^3; %y(k)=1.1*y(k-1)-0.2*sin(t(k-1))*y(k-1)^3; %y(k)=1.3*y(k-1)-0.6*sin(t(k-1))*y(k-1)^3; end plot(T,z,r) hold on plot(t,y,t,y,r*) title(Euler Method of Solving Initial Value Problem) legend(Integral Curve,Euler Curve) 附　应用程序　Euler_Method_File.m function varargout = odefile1(t,y,flag) switch flag case % Return dy/dt = f(t,y). varargout{1} = f(t,y); case init % Return default [tspan,y0,options]. [varargout{1:3}] = init; case jacobian % Return Jacobian matrix df/dy. varargout{1} = jacobian(t,y); otherwise error([Unknown flag flag .]); end % ------------------------------------------------- function dydt = f(t,y) dydt =[y(:)-2.*sin(t(:)).*y(:).^3]; % ------------------------------------------------- function [tspan,y0,options] = init tspan =[0;3