常微分方程数值解法

常微分方程数值解法

【作用】微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步：

1. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。

2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。

3. 运用这些规律列出方程和定解条件。基本模型

1. 发射卫星为什么用三级火箭

2. 人口模型

3. 战争模型

4. 放射性废料的处理通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来” 的于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。

1. 改进Euler 法：

2. 龙格—库塔( Runge—Kutta )方法：

【源程序】

1. 改进Euler 法：

function [x,y]=eulerpro(fun,x0,x1,y0,n);%fun 为函数，(xO, x1)为x 区间，yO 为初始值，n 为子

区间个数

if nargin<5,n=5O;end

h=(x1-xO)/n;

x(1)=xO;y(1)=yO;

for i=1:n

x(i+1)=x(i)+h;

y1=y(i)+h*feval(fun,x(i),y(i));

y2=y(i)+h*feval(fun,x(i+1),y1);

y(i+1)=(y1+y2)/2;

end

调用command 窗口

f=i nlin e('-2*y+2*x A2+2*x') [x,y]=eulerpro(f,O,,1,1O)

2

x +2x , (0 < x < , y(0) = 1

求解函数y'=-2y+2

2. 龙格—库塔( Runge—Kutta )方法：

[t,y]=solver('F',tspan ，y0)

这里solver为ode45, ode23, ode113，输入参数F是用M文件定义的微分方程y'= f (x, y)右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。

注：ode45和ode23变步长的，采用Runge-Kutta算法。

ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长(变步长)的常微分方程数值解法，其整体截断误差为(△ 口人5解

决的是Nonstiff(非刚性)常微分方程。

ode45是解决数值解问题的首选方法，若长时间没结果，应该就是刚性的，可换用ode23 试试例如：咛二y十

odefun=@(t,y) (y+3*t)/t A2; % 定义函数

tspa n=[1 4]; %求解区间

y0=-2; %初值

[t,y]=ode45(odefu n,tspa n,y 0);

plot(t,y) % 作图

title('tA2y”=y+3t,y(1)=-2,1

lege nd('tA2y''=y+3t')

xlabel('t')

ylabel('y')

%精确解

dsolve('tA2*Dy=y+3*t','y(1)=-2')

ans = (3*Ei(1,-1/t)+(-3*exp(-1)*Ei(1,-1)-2)/exp(-1))*exp(-1/t)

求解高阶常微分方程

需要求解的高阶常微分方程：

y r>= -iy + 3sin2f

求解的关键是将高阶转为一阶,odefun的书写.

F(y,y',y"...y(n-1),t)=0用变量替换,y1=y,y2=y'...注意odefun方程定义为行向量程序：

fun cti on Testode45

tspa n=[ ]; %求解区间

y0=[2 8]; % 初值

[t,x]=ode45(@odefu n,tspa n,y0);

plot(t,x(:,1),'-o',t,x(:,2),'-*')

lege nd('y1','y2')

title('y” ”=-t*y + eAt*y” +3si n2t')

xlabel('t')

ylabel('y')

end

fun cti on y=odefu n(t,x)

y=zeros(2,1); % 列向量

y(i)=x(2)；

y(2)=-t*x(1)+exp(t)*x (2)+3*si n(2*t);

end

【原理】

改进 Euler 法：

叫改进Euler 法。为了编程方便，将式子改写成

龙格一库塔（Runge — Kutta ）方法：

在区间［Xn, Xn+1］内多取几个点，将它们的斜率加权平均作为 K ,就有可能构造

出精度更高的计算公式。这就是龙格一库塔方法的基本思想。

要进-步捉高椿度*必绩取更多的点.如联4点构造如下形式的公式；

］儿+1 —儿+肌人杠+心鸟卜久怡+為孔）

*1 =于（兀”几）

■ H 一/（%亠◎再儿+庆叭'

（⑺

若=f （x a +a 2h,儿 +0』化 +00 札）

忑=/（%十吗乩儿十/V 叭十民哄十0押J

其中待定条数人G 厂久共13个，经过与推导2阶龙格一库塔公式类似、但也更杂的 计算.得到使局部误差y 仗曲）-y* =。（护）的11个方程•収既满足这些方程、又较 简单的一齟心勺，冋‘可得卜=几€尼小"八"几川軌式称为由

Euler 公式和梯形公式得到的预测一校正系统,