第七章 ODE 的数值解

1.使用ODE23 和ODE45 求解一阶方程

（1）要求解 ODE 的数值解，首先得定义表示方程的的函数。

（2）创建一个.m 文件，输入下面的内容

function ydot = eq1(t,y)

ydot = cos(t);

（3）通过调用 ODE32 函数来求解 ODE。这个函数使用使用二阶和三阶 Runge-Kutta（龙格-库塔）法对微分方程进行积分。

[t,y] = ode23(‘func_name’, [start_time, end_time], y(0));

[t,y] = ode23(‘eq1’,[0 2*pi],2);

f = 2 + sin(t);

plot(t,y,‘o’,t,f),xlabel(‘t’),ylabel(‘y(t)’),axis([0 2*pi 0 4]);//绘图

（4）把误差函数起名为 err。首先我们为储存误差点的数组分配内存空间，我们使用一个所有元素都为零的数组来表示：

err = zeros(size(y));//配置误差函数

（5）for 循环遍历数据，计算每个点上的相对误差：

for i = 1:1:size(y)

err(i) = abs((f(i)-y(i))/f(i));//相对误差计算

end

err

err =

1.0e-003 *

0

0.0001

0.0091

0.0301

0.0743

0.2115

0.2876

0.3780

0.4659

0.5597

0.6480

0.6907

0.6050

0.4533

0.3164

0.2414

0.2129

emax = max(err)//获取最大误差

emax =6.9075e-004

ode45 函数使用更高阶 Runge-Kutta 公式。

[t,w] = ode45(‘eq1’,[0 2*pi],2);//基于之前脚本文件eq1

f = 2 + sin(t);

plot(t,w,‘o’,t,f),xlabel(‘t’),ylabel(‘y(t)’),axis([0 2*pi 0 4])//绘图

比较解的大小：

size(w),size(y)

ans = 45 1

ans = 17 1

err = zeros(size(w));误差函数配置

for i = 1:1:size(w)

err(i) = abs((f(i)-w(i))/f(i));

end

wmax = max(err)//最大误差

wmax = 4.9182e-006

emax/wmax

ans =140.4473

结论：ode32比 ode45 得到的相对误差大得多，实际上它差不多大 141 倍。如果误差对你非常重要，那么使用 ode45 将是明智的选择。

2.二阶方程求解

（1）配置脚本

function xdot = eqx(t,x);

%创建数组储存数据

xdot = zeros(2,1);

xdot(1) = -x(1)^2 + x(2);

xdot(2) = -x(1) - x(1)* x(2);

（2）求解函数

[t,x] = ode45(‘eqx’,[0 10],[0,1]);

（3）绘图

plot(t,x(:,1),t,x(:,2),‘–’),xlabel(‘t’), axis([0 10 -1.12 1.12])

（4）显示两个函数的图象

plot(x(:,1),x(:,2))