数学实验“欧拉数值算法,rungekutta数值算法”实验报告(内含matlab程序).doc

西京学院数学软件实验任务书 课程名称 数学软件实验 班级 数 0901 学号 0912020107 姓名 李亚强 实验课题 欧拉数值算法(显式，隐式，欧拉预估-校正法) ， Runge-Kutta 数值算法 实验目的 熟悉欧拉数值算法(显式，隐式，欧拉预估-校正法) ， Runge-Kutta 数值算法 实验要求 运用Matlab/C/C++/Java/Maple/Mathematica 等其中 一种语言完成 实验内容 欧拉数值算法(显式，隐式，欧拉预估-校正法) Runge-Kutta 数值算法 成绩 教师- 1 - 实验二十四实验报告 1、实验名称：欧拉数值算法(显式，隐式，欧拉预估-校正法) ， Runge-Kutta 数值算法。 2、实验目的：进一步熟悉欧拉数值算法(显式，隐式，欧拉预估- 校正法) ，Runge-Kutta 数值算法。 3、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica 等其 中一种语言完成程序设计。 4、实验原理： 1．欧拉数值算法(显式)： 微分方程里最简单的方程形式莫过于一阶常微分方程的 初值问题，即： 0 ( , ) ( ) dy f x y a x b dx y a y          其中 ， 为常数。 a b 因为其简单但又是求解其他方程的基础，所以发展了许 多典型的解法。所有算法中的 就是代表上式中 ，而 f ( , ) f x y 表示 ， 表示 。 y f ( , ) f x y y   x f ( , ) f x y x   简单欧拉法是一种单步递推算法。简单欧拉法的公式如- 2 - 下所示： 1 ( , ) n n n n y y hf x y    简单欧拉法的算法过程介绍如下: 给出自变量 的定义域 ，初始值 及步长 。 x [ , ] a b 0 y h 对 ，计算 0,1 ,( ) / k b a h    1 ( , ) k k k k y y hf x y    2．欧拉数值算法(隐式)： 隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示： 1 1 1 ( , ) n n n n y y hf x y      隐式欧拉法是一阶精度的方法，比它精度高的公式是： 1 1 1 [ ( , ) ( , )] 2 n n n n n n h y y f x y f x y       隐式欧拉的算法过程介绍如下。 给出自变量 的定义域 ，初始值 及步长 。 x [ , ] a b 0 y h 对 ，用牛顿法或其他方法求解方 0,1 ,( ) / k b a h    1 1 1 ( , ) k k k k y y hf x y      得出 。 1 k y  3．欧拉预估-校正法： 改进欧拉法是一种二阶显式求解法，其计算公式如下所 示：- 3 - 1 [ , ( , )] 2 2 n n n n n n h h y y hf x y f x y      1 1 ( , ) [ ( , ) ( , )] 2 n n n n n n n n t y hf x y h y y f x y f x t             4．Runge-Kutta 数值算法： 二阶龙格-库塔法有多种形式，除了改进的欧拉法外， 还有中点法。中点法计算公式为： 1 [ , ( , )] 2 2 n n n n n n h h y y hf x y f x y      5、实验内容： %欧拉数值算法(显式) function [h,k,X,Y,P,REn]=Qeuler1(funfcn,x0,y0,b,n,tol) x=x0; h=(b-x)/n; X=zeros(n,1); y=y0; Y=zeros(n,1); k=1; X(k)=x; Y(k)=y ; for k=2:n+1fxy=f(funfcn,x,y); delta=norm(h*fxy, inf );wucha=tol*max(norm(y, inf ),1.0); if delta>=wucha x=x+h; y=y+h*fxy; X(k)=x;Y(k)=y ; end plot(X,Y, rp ) grid label( 自变 量 X ), ylabel( 因变量 Y ) title( 用向前欧拉(Euler)公式计算dy/dx=f(x,y) ，y(x0)=y0在 [x0,b]上的数值解 ) end P=[X,Y]; syms dy2, REn=0.5*dy2*h^2;- 4 - %欧拉数值算法(隐式) function [X,Y,n,P]=Heuler1(funfcn,x0,b,y0,h,tol) n=fix((b-x0)/h); X=zeros(n+1,1); Y=zeros(n+1,1); k=1; X(k)=x0; Y(k,:)=y0; Y1(k,:)=y0; clc,x0,h,y0 for i=2:n+1 X(i)=x0+h; Y(i,:)=y0+h*f(funfcn,x0,y0); Y1(i,:)=y0+h*f(funfcn,X(i),Y(i,:)); Wu=abs(Y1(i,:)-Y(i,:)); while Wu>tolp=Y1(i,:); Y1(i,:)=y0+h*f(funfcn,X(i),p); Y(i,:)=p;end x0=x0+h; y0=Y1(i,:); Y(i,:)=y0; plot(X,Y, ro ) grid on xlabel( 自变量 X ), ylabel( 因变量 Y ) title( 用向后欧拉公式 计算dy/dx=f(x,y)，y(x0)=y0 在[x0,b] 上 的数值解 ) endX=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1; P=[n ,X,Y] %欧拉预估-校正法 function [H,X,Y,k,h,P]=gaiEuler(funfcn,x0,b,y0,tol) pow=1/3; if nargin x)if x+h>b h=b-x;endfxy=f(funfcn,x,y); fxy=fxy(:); delta=norm(h*fxy, inf ); wucha=tol*max(norm(y, inf ),1.0); if deltalength(X) X=[X;zeros(p,1)]; Y=[Y;zeros(p,length(y))]; H=[H;zeros(p,1)];end H(k)=h;X(k)=x;Y(k,:)=y ; plot(X,Y, mh ), grid xlabel( 自变量 X ), ylabel( 因变量 Y ) title( 用改 进的欧拉公式计算dy/dx=f(x,y)，y(x0)=y0 在[x0,b] 上的数值解 ) end if delta~=0.0 h=min(h*8,0.9*h*(wucha/delta)^pow); end end if