西京学院数学软件实验任务书24

1、西京学院数学软件实验任务书课程名称数学软件实验班级数0901学号0912020119姓名王震实验课题欧拉数值算法（显式，隐式，欧拉预估-校正法），Runge-Kutta数值算法实验目的熟悉欧拉数值算法（显式，隐式，欧拉预估-校正法），Runge-Kutta数值算法实验要求运用Matlab/C/C+/Java/Maple/Mathematica等其中一种语言完成实验内容欧拉数值算法（显式，隐式，欧拉预估-校正法）Runge-Kutta数值算法成绩教师【实验课题】欧拉数值算法（显式，隐式，欧拉预估-校正法），Runge-Kutta数值算法【实验目的】熟悉欧拉数值算法（显式，隐式，欧拉预估-校正法）

2、，Runge-Kutta数值算法【实验内容】1、欧拉数值算法对于微分方程： （1）的数值解，当函数对于满足 Lipschitz 条件：时，初值问题（1）存在唯一解。数值解法就是寻求解在离散点上的近似值，并且有，将向前差商：代入(1)并记，由此可得显式欧拉公式： (2)若使用向后差商：和(1)，可得隐式欧拉公式： (3)对方程的两端从到积分得： (4)将梯形公式应用于上式右端的积分项，可得梯形公式： (5)上式可视为显式欧拉公式和隐式欧拉公式的算术平均。 实际计算中，可将显式欧拉公式和梯形公式结合使用，建立欧拉预估校正公式，或改进的欧拉公式：(6)即：2、Runge-Kutta数值算法为了进一步

3、提高精度，在上可取多个点，预报相应点的斜率值，对这些斜率值加权平均作为平均斜率值。利用泰勒展开，比较相应系数，从而确定在尽可能高的精度下有关参数应满足的条件。较常用的有三阶龙格-库塔公式：【程序】%向前欧拉数值算法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+1 fxy=feval(funfcn,x,y); delta=norm(h*fxy,'inf'

4、); wucha=tol*max(norm(y,'inf'),1.0);if delta>=wuchax=x+h; y=y+h*fxy; X(k)=x;Y(k)=y'endplot(X,Y,'rp')gridlabel('自变量 X'), ylabel('因变量 Y')title('用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解') endP=X,Y; syms dy2,REn=0.5*dy2*h2;%向后欧拉数值算法function X,Y,n,P=H

5、euler1(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;for i=2:n+1X(i)=x0+h; Y(i,:)=y0+h*feval(funfcn,x0,y0);Y1(i,:)=y0+h*feval(funfcn,X(i),Y(i,:); Wu=abs(Y1(i,:)-Y(i,:);while Wu>tol p=Y1(i,:); Y1(i,:)=y0+h*feval(funfcn,X(i),p); Y(i,:)=p;end

6、x0=x0+h; y0=Y1(i,:); Y(i,:)=y0; plot(X,Y,'ro')grid onxlabel('自变量 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 E=MendEuler(f,a,b,n,ya)% f:微分方程右端函数句柄% a,b:自变量取值区间的两个端点% n:区间等分的个数% y

7、a:函数初值y(a)% E=x',y':自变量X 和解Y 所组成的矩阵h=(b-a)/n;y=zeros(1,n+1);x=zeros(1,n+1);y(1)=y(a);x=a:h:b;for i=1:ny1=y(i)+h*feval(f,x(i),y(i);y2=y(i)+h*feval(f,x(i+1),y1);y(i+1)=(y1+y2)/2;endE=x',y'%三阶龙格库塔function k,X,Y,fxy,wch,wucha,P=RK3(funfcn,fun,x0,b,C,y0,h)x=x0; y=y0;p=128; n=fix(b-x0)/h);

8、fxy=zeros(p,1);wucha=zeros(p,1);wch=zeros(p,1); X=zeros(p,1); Y=zeros(p,length(y); k=1; X(k)=x; Y(k,:)=y'% 绘图.clc,x,h,y %计算 %fxy=fxy(:);for k=2:n+1x=x+h;a2=C(4); a3=C(5);b21=C(6); b31=C(7); b32=C(8);c1=C(1); c2=C(2); c3=C(3); x1=x+a2*h; x2=x+a3*h; k1=feval(funfcn,x,y); y1=y+b21*h*k1;k2=feval(fun

9、fcn,x1,y1); y2=y+b31*h*k1+b32*h*k2; k3=feval(funfcn,x2,y2);fxy(k)=feval(fun,x); y=y+h*(c1*k1+c2*k2+c3*k3); X(k)=x; Y(k,:)=y; k=k+1;plot(X,Y,'rp',X,fxy,'bo'), grid,xlabel('自变量 X'), ylabel('因变量 Y')legend('用三阶龙格-库塔方法计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解','y/dx=f(x,y)，y(x0)=y0的精确解y=f(x)')end%计算误差.for k=2:n+1wucha(k)=norm