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

y(1)=y0; h=(xn-x0)/n; x=x0:h:xn; disp('改进欧拉法结果如下：') disp(['y','(',num2str(x0),')','=',num2str(y0)]); for i=2:n+1 %利用循环语句实现欧拉法利用第一个节点的 函数值计算 k(1)=f(x0,y(i-1)); k(2)=f(x0,y(i-1)+h*k(1)); y(i)=y(i-1)+h*(1/2*k(1)+1/2*k(2)); %随后所有节点函数值的过 程. x0=x0+h; disp(['y','(',num2str(x0),')','=',num2str(y(i))]); end plot(x,y,'k*-'); %利用求出的x坐标和y坐标画出解的近 似图形. xlabel('x') ylabel('y') title('改进欧拉法求出的折线') hold on %-------若原微分方程的理论解能求出，则可画出积分曲线比较-------------x=0:0.01:0.5; y=x+exp(-x); plot(x,y,'b-') %-----------------------------------------------------------------legend('改进欧拉法求出的折线','解曲线') hold off function y=f(x,z) y=-z+x+1; 对于例1的结果为 输出图像为
实验七 常微分方程数值解法
实验目的
1、 掌握常微分方程数值解法中欧拉法和改进欧拉法的基本原理 和算法； 2、 培养Matlab数学软件的编程与上机调试能力.
实验课时
4课时
实验准备
初值问题 的解析解法在常微分方程中有介绍,但是实际中遇到的常微分方程一般 很难得到解析解,所以我们一般用数值方法求出它的数值解. 数值解法的 结果是自变量在一系列离散点处的函数的近似值,一般情况下取步长为 定数. 1.欧拉法 利用一阶微分公式 可得差分方程 其中, , ,则得到初值问题的数值解. 2.后退欧拉法 利用一阶微分公式 可得差分方程 此为后退欧拉法. 4. 梯形求积公式 利用数值积分 可得梯形公式 5. 改进欧拉公式 一般来讲,隐式格式要比显示格式具有较好的数值稳定性,所以常常被 采用,但是为了避免隐式格式的计算量,一般采用预估-校正的方法.对于梯 形公式,用欧拉公式预估再用梯形公式校正,可得改进欧拉公式:
来自百度文库 欧拉法和改进欧拉法的效果比较图为
实验习题
1. 利用欧拉法程序,求解初值问题 （注意：若理论解未知，解曲线画不出来，应删掉相关程序语句） (1) (2) 并将欧拉法求出的折线和改进欧拉法求出的折线放在同一个图中，比较 效果。 2. 某跳伞者在 时刻从飞机上跳出,假设初始时刻的垂直速度为0,且跳伞者垂直下落.已 知空气阻力为 ,其中c为常数,v为垂直速度,向下方向为正.可求得跳伞者的速度满足微分 方程方程 ,g为重力加速度.若此跳伞者质量为 ,且已知
,试试编程利用欧拉公式画出 的速度曲线(取 ).
实验算例
1. 欧拉法 利用matlab中的一个循环语句即可实现欧拉法中的利用第一个节点计
算随后所有节点的工作. function Euler %输出节点的x值和y值, clc %------输入(x0,xn)为求解区间,y0=y(x0)为初始条件,n为区间的 等分个数------x0=0; y0=1; xn=0.5; n=50; %-------------------------------------------------------------------y(1)=y0; h=(xn-x0)/n; x=x0:h:xn; disp('欧拉法结果如下：') disp(['y','(',num2str(x0),')','=',num2str(y0)]); for i=2:n+1 %利用循环语句实现欧拉法利用第一个节 点的函数值计算 y(i)=y(i-1)+h*f(x0,y(i-1)); %随后所有节点函数值的过程. x0=x0+h; disp(['y','(',num2str(x0),')','=',num2str(y(i))]); end plot(x,y,'ro-'); %利用求出的x坐标和y坐标画出解的 近似图形. xlabel('x') ylabel('y') title('欧拉法求出的折线') hold on %-------若原微分方程的理论解能求出，则可画出积分曲线比 较--------------x=0:0.01:0.5; y=x+exp(-x); plot(x,y,'b-') hold off %-----------------------------------------------------------------function y=f(x,z) y=-z+x+1;
例1. 用欧拉法求解初值问题, 并画出解图形. >> 输出图像为
另外,还可以改变输入的参数n,加密节点,得到更精确的数值解. 2. 改进欧拉法 function Gaijin_Euler %输出节点的x值和y值, clc close all %------输入(x0,xn)为求解区间,y0=y(x0)为初始条件,n为区间的等 分个数------x0=0; y0=1; xn=0.5; n=50; %--------------------------------------------------------------------