实验四常微分方程初值问题数值解法

时间：2025-07-12

实验报告

课程名称数值分析实验项目常微分方程问题初值问题数值解法

一. 实验目的

1、理解如何在计算机上实现用Euler法、改进Euler法、Runge－Kutta算法求一阶常微分方程初值问题

y (x) f(x,y),x [a,b]

y(a) y1

的数值解。

2、利用图形直观分析近似解和准确解之间的误差。 3、学会Matlab提供的ode45函数求解微分方程初值问题。

二、实验要求

（1）按照题目要求完成实验内容；（2）写出相应的Matlab 程序；

（3）给出实验结果(可以用表格展示实验结果)；（4）分析和讨论实验结果并提出可能的优化实验。（5）写出实验报告。三、实验步骤

1、用编好的Euler法、改进Euler法计算书本P167 的例1、P171例题3。（1）取h 0.1，求解初值问题

2x y(x) y ,x [0,1],

y(0) 1

（2）取h 0.1，求解初值问题

y (x) y x 1,x [0,0.5],

y(0) 1

2、用Runge－Kutta算法计算P178例题、P285实验任务（2）

（1）取h 0.1，求解初值问题

y (x) y2,x [0,0.5],

y(0) 1

（2）求初值问题

1 2

y(x) ( y x 4x 1),x [0,0.5],

y(0) 0

的解y(x)在xi ih(h 0.05)处的近似值yi，并与问题的解析解y(x) e相比较。

x2 1

3、用Matlab绘图函数plot(x,y)绘制P285实验任务（2）的精确解和近似解的图形。

4、使用matlab中的ode45求解P285实验任务（2），并绘图。

四、实验结果

1、Euler算法程序、改进Euler算法程序；

Euler算法程序：function [x,y]=euler_f(ydot_fun, x0, y0, h, N) % Euler（向前）公式，其中

% ydot_fun --- 一阶微分方程的函数 % x0, y0 --- 初始条件 % h --- 区间步长 % N --- 区间的个数 % x --- Xn 构成的向量 % y --- Yn 构成的向量

x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:N

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

y(n+1)=y(n)+h*feval(ydot_fun, x(n), y(n)); end

改进Euler算法程序：

function [x,y]=euler_r(ydot_fun, x0, y0, h, N) % 改进Euler公式，其中

% ydot_fun --- 一阶微分方程的函数 % x0, y0 --- 初始条件 % h --- 区间步长 % N --- 区间的个数 % x --- Xn 构成的向量 % y --- Yn 构成的向量

x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:N

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

ybar=y(n)+h*feval(ydot_fun, x(n), y(n));

y(n+1)=y(n)+h/2*(feval(ydot_fun, x(n), y(n))+feval(ydot_fun, x(n+1), ybar)); end

2、用Euler算法程序、改进Euler算法求解P167例题1的运行结果；（1.）Euler算法程序： x =

Columns 1 through 8

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.7000

Columns 9 through 11

0.8000 0.9000 1.0000 y =

0.6000

Columns 1 through 8

1.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5090 1.5803

Columns 9 through 11

1.6498 1.7178 1.7848

（2）改进Euler算法： x =

Columns 1 through 7

0 0.1000 0.2000 0.3000 0.6000

Columns 8 through 11

0.7000 0.8000 0.9000 1.0000 y =

Columns 1 through 7

1.0000 1.0959 1.1841 1.2662 1.4860

0.4000 0.5000 1.3434 1.4164

Columns 8 through 11

1.5525 1.6165 1.6782 1.7379

3、Runge－Kutta算法程序；

function [x,y]=Runge_Kutta4(ydot_fun, x0, y0, h, N) % 标准四阶Runge_Kutta公式，其中 % ydot_fun --- 一阶微分方程的函数 % x0, y0 --- 初始条件 % h --- 区间步长 % N --- 区间的个数 % x --- Xn 构成的向量 % y --- Yn 构成的向量

x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; for n=1:N

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

k1=h*feval(ydot_fun, x(n), y(n));

k2=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k1); k3=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k2); k4=h*feval(ydot_fun, x(n)+h, y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end

4、用Runge－Kutta算法求解P178例题、P285实验任务（2），计算结果如下（其中yi表示数值解，y(xi)表示解析解，结果保留八位有效数字）：求解P178例题： x =

0 0.1000 0.2000 0.3000 0.4000 0.5000 y =

1.0000 1.1111 1.2500 1.4286 1.6667 2.0000

5、P285实验任务（2）精确解与近似解的图形比较

6、用matlab中的ode45求解P285实验任务（2） ans =

Columns 1 through 7

0 0.0001 0.0002 0.0003 0.0004 0.0009 0.0014 0 -0.0001 -0.0001 -0.0002 -0.0002 -0.0005 -0.0007

Columns 8 through 14

0.0019 0.0024 0.0049 0.0074 0.0099 0.0125 0.0250 -0.0010 -0.0012 -0.0024 -0.0037 -0.0049 -0.0061 -0.0118

Columns 15 through 21

0.0375 0.0500 0.0625 0.0750 0.0875 0.1000 0.1125 -0.0172 -0.0222 -0.0268 -0.0312 -0.0351 -0.0388 -0.0420

Columns 22 through 28