实验七、常微分方程数值解法

实验七常微分方程数值解法

实验目的

1.掌握常微分方程数值解法的基本原理和算法. 2.通过实验掌握各个算法的MATLAB程序.

实验课时

4课时

实验准备

初值问题

?dy??f?x,y?,a?x?b ?dx?y?a????的解析解法在常微分方程中有介绍,但是实际中遇到的常微分方程一般很难得到解析解,所以我们一般用数值方法求出它的数值解. 数值解法的结果是自变量x在一系列离散点

x1?x2??,xn,?处的函数f?x?的近似值y1,y2,?,yn,?,一般情况下取步长h?xn?1?xn为定数.

1.欧拉法

利用一阶微分公式

y?xn?1??y?xn?dy ?dxx?xnh可得差分方程

yn?1?yn?hf?xn,yn?,n?0,1,...,N?1

其中,y0??,yn?y?xn?,则得到初值问题的数值解.

2.后退欧拉法

利用一阶微分公式

dydx可得差分方程

?x?xn?1y?xn?1??y?xn?

h??yn?1?yn?hf?xn?1,yn?1?,n?0,1,...,N?1 ?y????0此为后退欧拉法.

4. 梯形法

利用数值积分

y?xn?1??y?xn???可得梯形公式

xn?1xnf?x,y?dx?h?f?xn,y?xn???f?xn?1,y?xn?1????? 2h??yn?1?yn???f?xn,yn??f?xn?1,yn?1???,n?0,1,...,N?1

2???y0??5. 改进欧拉法

一般来讲,隐式格式要比显示格式具有较好的数值稳定性,所以常常被采用,但是为了避免隐式格式的计算量,一般采用预估-校正的方法.对于梯形公式,用欧拉公式预估再用梯形公式校正,可得改进欧拉公式:

实验算例

1. 欧拉法

利用matlab中的一个循环语句即可实现欧拉法中的利用第一个节点计算随后所有节点的工作. function [x,y]=Euler(y0,a,b,n) clc y(1)=y0; h=(b-a)/n; x=a:h:b; for i=2:n+1 y(i)=y(i-1)+h*f(x(i-1),y(i-1)); end plot(x,y,'*'); 例1. 用欧拉法求解初值问题, 并画出解图形.

????

=???+??+1 ????,0≤??≤0.5 ?? 0 =1

首先存储 function z=f(x,y) z=-y+x+1; 命令窗口输入 >>[x,y]=Euler(1,0,0.5,10) 可得输出图像：

另外,还可以改变输入的参数n,加密节点,得到更精确的数值解。若能求出精确解，则可以画图观察精确程度。

放大之后看效果

2. 改进欧拉法 function [x,y]=Gaijin_Euler(y0,a,b,n) clc y(1)=y0; h=(b-a)/n; x=a:h:b; for i=2:n+1 k(1)=f(x(i-1),y(i-1)); k(2)=f(x(i-1),y(i-1)+h*k(1)); y(i)=y(i-1)+h*(1/2*k(1)+1/2*k(2)); end plot(x,y,'d'); 对于例1的结果为

两个方法画图比较

增加节点效果

实验习题

1. 编写梯形法程序，注意梯形法是隐式格式，需要用迭代法求解 ????+1=????+???(????+????)

(??+1) ?(??)

????+1=????+[?? ????+???? +??(????+????+1)]

22. 利用欧拉法，改进欧拉法，梯形法程序,求解初值问题

????

=?????? ????,0≤??≤1 ?? 0 =?1

并将三个方法放在同一个图中，和精确解比较。

(0)