《微分方程在科技》word版.doc

高等教育出版社 教育电子音像出版社 作者：任玉杰 第十章 常微分方程（组）求解第十章 常微分方程（组）求解微分方程在科技、工程、经济管理以及生态、环境、人口、交通等各个领域中常用于建立数学模型，它是研究函数变化规律的有力工具.如在研究弹性物体的振动，电阻、电容、电感电路的瞬变，热量在介质中的传播，抛射体的轨迹，以及污染物浓度的变化，人口增长的预测，种群数量的演变，交通流量的控制等等过程中，作为研究对象的函数，要和函数的导数一起，用一个符合其内在规律的方程，即微分方程来描述.建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验.如果能得到解析解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解，而绝大多数变系数方程、非线性方程都很难得到解析解.另外，有时即使能求出解析解，也会由于很难从解析解中计算函数的值而不实用.例如，容易求出初值问题的解为eed.但是，对于给定的，要计算函数值还需要用数值积分的方法.于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段.为了将微分方程的数值解与解析解比较，本章首先简单介绍用MATLAB解微分方程的解析解（符号解）的方法，然后重点介绍用MATLAB程序计算微分方程的数值解法.10.1 常微分方程（组）的MATLAB符号求解在MATLAB系统中提供了常微分方程（组）符号解（Symbolic solution of ordinary differential equations）的函数dsolve.在调用此函数前，必须首先将给定的常微分方程（组）中的一阶导数用D表示，例如，写成Dy.阶导数用Dn表示，例如，表示为D6y.由此，常微分方程写成DD.下面具体介绍函数dsolve解常微分方程（组）的方法.10.1.1 用MATLAB求常微分方程（组）的通解用MATLAB函数dsolve求常微分方程 （10.1）的通解的主要调用格式如下：调用格式一： S=dsolve (eqn，var)输入的量：eqn是常微分方程（10.1）改用符号方程表示的常微分方程DDD.Var表示自变量，默认的自变量为t.输出的量：S是常微分方程（10.1）的通解.求常微分方程组 （10.2）的通解的主要调用格式如下：调用格式二： S=dsolve (eqn1,eqn2, . ,eqnm，var)输入的量：eqn1, eqn2 , . , eqnm 分别是常微分方程组（10.2）中用符号方程表示的个常微分方程.默认的自变量为t，也可以将自变量t变为其它的符号变量var.输出的量：S是常微分方程组（10.2）的通解.10.1.2 用MATLAB求常微分方程（组）的特解如果给定常微分方程（10.1）的初始条件 （10.3）则求方程（10.1）的特解的主要调用格式如下：调用格式三： S=dsolve (eqn，condition1，conditionn，var)输入的量：eqn是常微分方程（10.1）改用符号方程表示的常微分方程DDD；condition1，conditionn是初始条件（10.3）；Var表示自变量，默认的自变量为t.输出的量：S是常微分方程（10.1）的特解.同理，如果给定常微分方程组（10.2）的初始条件，则求方程（10.2）的特解的主要调用格式如下： 调用格式四：S=dsolve(eqn1,eqn2, . ,eqnm，condition1，condition2，var)输入的量：eqn1, eqn2, . ,eqnm 分别是常微分方程组（10.2）中用符号方程表示的个常微分方程；condition1，condition2，是常微分方程组（10.2）的初始条件；默认的自变量为t，也可以将自变量t变为其它的符号变量var.输出的量：S是常微分方程组（10.2）的特解.10.1.3 线性常微分方程组的解法如果线性常微分方程组中的未知量较多，用上面的求解方法就不方便，这时我们可以用下面介绍的求解齐次线性常微分方程组和非齐次线性常微分方程组的方法求解.求解齐次线性常微分方程组的MATLAB主程序输入量：A是齐次线性常微分方程组的矩阵.输出量：X是的解，E是矩阵的特征值，V是矩阵的特征向量.名为qcxxcwz.m的求解齐次线性常微分方程组的MATLAB主程序如下：function X,E,V=qcxxcwz(A)syms xE=eig(A);V,n=eig(A); X=exp(E*x)*V;10.3 欧拉（Euler）方法的MATLAB程序欧拉方法是一种最简单的解微分方程的数值方法，其基本想法与数值微分法基本思想相同，是在小区间上用差商代替（10.5）中方程左端的导数，而方程右端的已知函数中的在小区间上取值不同，则有以下不同的方法.本节主要介绍用向前欧拉公式计算初值问题的三种MATLAB程序和具体的算例.10.3.1 向前欧拉公式及其误差估计如果中的取小区间的左端点，由上节用三种不同的方法得到相同的向前欧拉公式 （10.10）欧拉（Euler）方法就是用差分方程初值问题的解作为常微分方程初值问题（10.5）的数值解，即将初值代入式（10.7），逐次求出.例10.2.1 用欧拉方法求初值问题的数值解，分别取，并计算误差，画出精确解和数值解的图形.解 编写并保存名为Eulerli1.m的MATLAB计算和画图的主程序如下function P=Eulerli1(x0,y0,b,h)n=(b-x0)/h; X=zeros(n,1); Y=zeros(n,1); k=1; X(k)=x0; Y(k)=y0; for k=1:nX(k+1)=X(k)+h; Y(k+1)=Y(k)+h*(X(k)-Y(k); k=k+1;endy=X-1+2*exp(-X); plot(X,Y,mp,X,y,b-)gridxlabel(自变量 X), ylabel(因变量 Y)title(用向前欧拉公式求dy/dx=x-y，y(0)=1在0,1上的数值解和精确解y=x-1+2 exp(-x) legend(h=0.075时，dy/dx=x-y，y(0)=1在0,1上的数值解,精确解y=x-1+2 exp(-x)jwY=y-Y;xwY=jwY./y;k1=1:n;k=0,k1; P=k,X,Y,y,jwY,xwY;在MATLAB工作窗口输入下面的程序x0=0;y0=1;b=1;h=0.0750;P=Eulerli1(x0,y0,b,h)运行后屏幕显示取h=0.075 0时，该初值问题的精确解y、用欧拉方法求此初值问题与自变量X对应的数值解Y，Y的绝对误差jwY和相对误差xwY（略），精确解y和数值解Y的图形（见图102）.图102 取h=0.075 0时的精确解和数值解 图103 取h=0.007 5时的精确解和数值解在MATLAB工作窗口输入下面的程序h1=0.0075; P1=Eulerli1(x0,y0,b,h1)legend(h1=0.0075时，dy/dx=x-y，y(0)=1在0,1上的数值解,精确解y=x-1+2 exp(-x)运行后屏幕显示取h=0.007 5时，该初值问题的精确解y、用欧拉方法求此初值问题与自变量X对应的数值解Y，Y的绝对误差jwY和相对误差xwY（略），精确解y和数值解Y的图形（见图103）.由图102和图103可以看出，h的绝对值越小，数值解Y与精确解y的误差越小.10.3.2 向前欧拉方法的三种MATLAB程序 下面分别介绍三种用向前欧拉公式求解常微分方程初值问题（10.5）的数值解及其截断误差公式的MATLAB程序，并且通过例题说明这三种程序的应用.(一) 向前欧拉公式的MATLAB主程序1下面介绍根据向前欧拉公式（10.8），现提供的数值计算常微分方程初值问题（10.5）的MATLAB程序.用向前欧拉公式（10.8）求解常微分方程初值问题（10.5）的数值解及其截断误差公式的MATLAB主程序1输入量：funfcn是函数，x0和y0是初值b是自变量的最大值，n是自变量取值的个数，tol是精度.输出量：数组 X的元素是由自变量的取值组成，数组Y的元素是利用向前欧拉公式（10.8）求出常微分方程初值问题（10.5）在X的元素处的数值解，REn是在每个子区间上的局部截断误差公式，h是步长，其中dy2是的2阶导数.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); 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;例10.3.2 用向前欧拉公式（10.8）求解初值问题，分别取，并将计算结果与精确解作比较，写出在每个子区间上的局部截断误差公式，画出数值解与精确解在区间上的图形.图104 =10,100时，在0,1上的数值解和精确解解 输入程序 subplot(2,1,1) x0=0;y0=1;b=1-1.e-4;n=100;tol=1.e-4;h1,k1,x1,Y1,P1,Ren1=QEuler1(funfcn,x0,y0,b,n,tol) hold onS1= 8/3*x1-29/9+38/9*exp(-3*x1), plot(x1,S1,b-)title(用向前欧拉公式计算dy/dx=8x-3y-7，y(0)=1在0,1上的数值解)legend(n=100时，dy/dx=8x-3y-7，y(0)=1在0,1上的数值解, dy/dx=8x-3y-7，y(0)=1在0,1上的精确解)hold offjdwuc1=S1-Y1; jwY1=S1-Y1;xwY1=jwY1./S1;k1=1:n;k=0,k1;P1=k,x1,Y1,S1,jwY1,xwY1subplot(2,1,2)n1=10; h2,k2,x2,Y2,P2,Ren2=QEuler1(funfcn,x0,y0,b,n1,tol) hold onS1 = 8/3*x2-29/9+38/9*exp(-3*x2), plot(x2,S1,b-)legend(n=10时，dy/dx=8x-3y-7，y(0)=1在0,1上的数值解, dy/dx=8x-3y-7，y(0)=1在0,1上的精确解)hold offjwY2=S1-Y2;xwY2=jwY2./S1;k1=1:n1;k=0,k1; P2=k,x2,Y2,S1,jwY2,xwY2 运行后屏幕显示分别取时，所给的初值问题在上的自变量的值构成的数组Xi (i=1,2)，利用向前欧拉公式（10.8）求出的与Xi (i=1,2)对应的数值解Yi (i=1, 2)，Yi的绝对误差jwYi (i=1,2)和相对误差xwYi (i=1,2)（略），每个子区间上的局部截断误差公式Ren2 =1/200*dy2， 近似解Y与精确解的图形（见图104）. 由图104可以看出，取时，利用向前欧拉公式（10.8）求出的该初值问题在0,1上的数值解Yi与精确解y的绝对误差向量jwYi和相对误差向量xwYi较大，但是当=100时，数值解Yi的图形和精确解y的图形几乎重合，即节点的个数越多，误差越小.这说明当节点的个数足够大时，这种方法是十分有效的. (二) 向前欧拉公式的MATLAB主程序2下面介绍根据向前欧拉公式（10.8），现提供的数值计算常微分方程初值问题（10.5）的另外一种MATLAB程序.用向前欧拉公式（10.8）求解常微分方程初值问题（10.5）的数值解及其截断误差公式的MATLAB主程序2输入量：funfcn是函数，x0和y0是初值，b是自变量的最大值，h是步长.输出量：k是自变量取值的个数，向量X的元素是自变量的取值，向量Y的元素是利用向前欧拉公式（10.8）求出常微分方程初值问题（10.5）在向量X的元素处的数值解， REn是在每个子区间上的局部截断误差公式，其中dy2是的2阶导数.function k,X,Y,P,REn=Qeuler2(funfcn,x0,y0,b,h)x=x0; n=fix(b-x)/h); X=zeros(n+1,1); y=y0;Y=zeros(n+1,1); k=1; X(k)=x; Y(k)=y; for k=2:n+1X(k)=x+(k-1)*h,fxy=feval(funfcn,x,y), Y(k)=y+h*fxyy=Y(k), k=k+1,plot(X,Y,rp), grid,xlabel(自变量 X), ylabel(因变量 Y)title(用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解) endk1=1:n;k=0,k1; P=k,X,Y; syms dy2,REn=0.5*dy2*h2;例10.3.4 用向前欧拉公式（10.8）求解初值问题，取，并将计算结果与精确解作比较，写出在每个子区间上的局部截断误差公式，画出数值解与精确解在区间上的图形.解 输入程序 subplot(2,1,1) x0=0;y0=1;b=2;n=10;h=2/10;k,X,Y,P,REn=Qeuler2(funfcn,x0,y0,b,h)hold onS1=1/6*(6+12*X+30*exp(2*X).(1/2); plot(X,S1,b-)title(用向前欧拉公式计算dy/dx=y-x/(3y)，y(0)=1在0,2上的数值解)legend(n=10时，dy/dx=y-x/(3y)，y(0)=1在0,2上的数值解, dy/dx =y-x/(3y)，y(0)=1在0,2上的精确解)hold offjdwucY=S1-Y; jwY=S1-Y;xwY=jwY./Y;k1=1:n;k=0,k1;P1=k,X,Y,S1,jwY,xwYsubplot(2,1,2)n1=100;h1=2/100;k,X1,Y1,P1,Ren1=Qeuler2(funfcn,x0,y0,b,h1) hold onS2=1/6*(6+12*X1+30*exp(2*X1).(1/2); plot(X1,S2,b-)legend(n=100时，dy/dx=y-x/(3y)，y(0)=1在0,2上的数值解, dy/dx =y-x/(3y)，y(0)=1在0,2上的精确解)hold offjwY1=S2-Y1;xwY1=jwY1./Y1;k1=1:n1;k=0,k1;P2=k,X1,Y1,S2,jwY1,xwY1运行后屏幕显示取时，此初值问题在上的自变量的向量X，X1，利用向前欧拉公式（10.8）求出的与X，X1对应的数值解Y，Y1及其绝对误差向量jwY，jwY1和相对误差向量xwY，xwY1（略），每个子区间上的局部截断误差公式Ren2 =1/200*dy2， 数值解Y与精确解的图形（见图107）.图107 取时，用向前欧拉公式求初值问题数值解和精确解的图形(三) 自适应向前欧拉公式的MATLAB主程序根据自适应方法和向前欧拉公式（10.8），现提供的数值计算常微分方程初值问题（10.5）的另一个MATLAB程序.用自适应向前欧拉公式求解常微分方程初值问题（10.5）的数值解的MATLAB主程序2输入量：funfcn是函数，x0和y0是初值b是自变量的最大值， tol是精度.输出量：向量X的元素是自变量的取值，向量Y的元素是利用向前欧拉公式（10.8）和自适应方法求出常微分方程初值问题（10.5）在向量X的元素处的数值解，向量H是步长h序列，n是自变量取值的序号，k是自变量的取值.并画出数值解向量Y的图形.function H,X,Y,k,h,P=QEuler(funfcn,x0,b,y0,tol)%初始化.pow=1/3; if nargin 5 | isempty(tol), tol = 1.e-6; end;if nargin 6 | isempty(trace), trace = 0; end;x=x0; h=0.0078125*(b-x); y=y0(:);p=128; H=zeros(p,1);X=zeros(p,1); Y=zeros(p,length(y); k=1; X(k)=x; Y(k,:)=y;% 绘图.clc,x,h,y end% 主循环.while (xx) if x+hbh=b-x; end %计算斜率. fxy=feval(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); endH(k)=h;X(k)=x;Y(k,:)=y; plot(X,Y,rp), gridxlabel(自变量 X), ylabel(因变量 Y)title(用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解)end%更新步长.if delta=0.0h=min(h*8,0.9*h*(wucha/delta)pow);endendif (xtol p=Y1(i,:); Y1(i,:)=y0+h*feval(funfcn,X(i),p); Y(i,:)=p; endx0=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上的数值解)end X=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1; P=n,X,Y例10.3.6 用向后欧拉公式求解区间上的初值问题 的数值解，取步长，并与精确解作比较，在同一个坐标系中做出图形.然后再取，观察数值解与精确解误差的变化，说明与误差的关系.解 输入程序S1=dsolve(Dy=8*x-3*y-7,y(0)=1,x) x0=0;y0=1; b=2;tol=1.e-1;subplot(2,1,1)h1=0.01; X1,Y1,n,P1=Heuler1(funfcn,x0,b,y0,h1,tol) hold onS2= 8/3*X1-29/9+38/9*exp(-3*X1), plot(X1,S2,b-)legend(h=0.01用向后欧拉公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解,dy/dx=8x-3y-7，y(0)=1在0,2上的精确解)hold offjuwY1=S2-Y1; xiwY1=juwY1./Y1; L=P1,S2,juwY1,xiwY1subplot(2,1,2)h=0.05; X,Y,n,P=Heuler1(funfcn,x0,b,y0,h,tol) hold onS1 = 8/3*X-29/9+38/9*exp(-3*X), plot(X,S1,b-)legend(h=0.05用向后欧拉公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解,dy/dx=8x-3y-7，y(0)=1在0,2上的精确解)hold offjuwY=S1-Y; xiwY=juwY./Y; L=P,S1,juwY,xiwY运行后屏幕显示用向后欧拉公式计算此初值问题在0,2上的自变量X处数值解Y和精确解S1及其图形（见图1010），步长H, Y的相对误差xiwY和绝对误差juwY(略) .图1010 取h=0.05时，用向后欧拉公式求初值问题的数值解和精确解的图形10.4 改进的欧拉方法的MATLAB程序下面主要介绍梯形公式、改进的欧拉方法和它们的误差估计得MATLAB程序，并通过例题说明它们的具体应用.10.4.2 梯形公式的两种MATLAB程序根据（10.15）式和（10.22）式现提供求解常微分方程初值问题（10.5）的MATLAB程序1和自适应MATLAB程序2.(一) 梯形公式的MATLAB程序用梯形公式求解常微分方程初值问题（10.5）的数值解的MATLAB主程序1输入量：funfcn是函数，x0和y0是初值b是自变量的最大值，h是步长， tol是精度.输出量：向量X的元素是自变量的取值，向量Y的元素是利用（10.15）式、（10.22）式和自适应方法求初值问题（10.5）在向量X的元素处的数值解，向量H是步长h序列，n是自变量取值的序号，k是自变量取值的序号，并画出数值解向量Y的图形.function X,Y,n,P= odtixing1(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; fx0y0=feval(funfcn,x0,y0); Y(i,:)=y0+h*fx0y0;fxiyi=feval(funfcn,X(i),Y(i,:); Y1(i,:)=y0+h*(fxiyi+fx0y0)/2; % 主循环.Wu=abs(Y1(i,:)-Y(i,:);while Wutol p=Y1(i,:), fxip=feval(funfcn,X(i),p);Y1(i,:)=y0+h*(fx0y0+fxip)/2,P1=Y1(i,:), Y(i,:)=p1; endx0=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上的数值解)end X=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1; P=n,X,Y例10.4.1 用梯形公式求解区间上的初值问题，取步长，精度为，并与精确解作比较，在同一个坐标系中画出图形.解 输入程序 x0=0;y0=1; b=2; tol=0.1; h=0.05;X,Yt,n,Pt=odtixing1(funfcn,x0,b,y0,h,tol)hold onS1=8/3*X-29/9+38/9*exp(-3*X); plot(X,S1,b-), hold offlegend(h=0.05,用梯形公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解,dy/dx=8x-3y-7，y(0)=1在0,2上的精确解)juwYt=S1-Yt; xiwYt=juwYt./Yt; Lt=Pt,S1,juwYt,xiwYt运行后屏幕显示取精度为，分别用梯形公式和向前欧拉公式求解此初值问题在区间上的自变量X处数值解Yi(i=t,q)和精确解S1，步长H, Yi的相对误差xiwYi和绝对误差juwYi (略) 及其数值解和精确解的图形（见图1011）.图1011用梯形公式求解区间上的初值问题 (二) 自适应梯形公式的MATLAB程序现提供求解常微分方程初值问题（10.5）的自适应MATLAB程序2如下：用自适应梯形公式求解常微分方程初值问题（10.5）的数值解的MATLAB主程序输入量：funfcn是函数，x0和y0是初值，b是自变量的最大值， tol是精度.输出量：向量X的元素是自变量的取值，向量Y的元素是利用（10.15）式、（10.22）式和自适应方法求出初值问题（10.5）在向量X的元素处的数值解，向量H是步长h的序列，n是自变量取值的序号，k是自变量取值的.并画出数值解向量Y的图形.function H,X,Y,k,h,P=odtixing2(funfcn,x0,b,y0,tol)%初始化.pow=1/3; if nargin 5 | isempty(tol), tol = 1.e-6; end;if nargin 6 | isempty(trace), trace = 0; end;x=x0; h=0.0078125*(b-x); y=y0(:);p=128; n=fix(b-x0)/h); H=zeros(p,1);X=zeros(p,1); Y=zeros(p,length(y); k=1; X(k)=x; Y(k,:)=y;% 绘图.clc,x,h,y % 主循环.while (xx) if x+hbh=b-x; end %计算斜率. fxy=feval(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); endH(k)=h;X(k)=x;Y(k,:)=y; plot(X,Y,go), gridxlabel(自变量 X), ylabel(因变量 Y)title(用自适应梯形公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解)end%更新步长.if delta=0.0h=min(h*8,0.9*h*(wucha/delta)pow);endendif (x x0=0;y0=1; b=2; tol=1.e-1;Ht,X,Yt,k,h,Pt= odtixing2(funfcn,x0,b,y0,tol),hold onS1 = 8/3*X-29/9+38/9*exp(-3*X), plot(X,S1,b-), hold offhold on,Hq,X,Yq,k,h,Pq=QEuler(funfcn,x0,b,y0,tol), hold offgridlegend(用自适应梯形公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解,dy/dx= 8x-3y-7， y(0)=1在0,2上的精确解,用向前欧拉公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解)title(用自适应梯形公式计算dy/dx=8x-3y-7,y(0)=1在0,2上的数值解)运行后屏幕显示取精度为，分别用梯形公式和向前欧拉公式求解此初值问题在区间上的自变量X处数值解Yi(i=t,q)和精确解S1及其图形（见图1012），步长H等 (略).图1012 自适应梯形公式和向前欧拉公式求解初值问题由图可见，用自适应梯形公式求解此初值问题的数值解的图形（绿色的小圆圈）几乎与精确解（蓝色的虚线）完全重合在一起，而用向前欧拉公式计算的数值解的图形（红色的小五角星）与精确解（蓝色的虚线）在比梯形公式的数值解误差大，所以，梯形公式比向前欧拉公式精度高.10.4.4 改进的欧拉公式的MATLAB程序现提供求解常微分方程初值问题（10.5）的MATLAB程序.用自适应改进的欧拉公式求解常微分方程初值问题（10.5）的数值解的MATLAB主程序输入量：funfcn是函数，x0和y0是初值b是自变量的最大值， tol是精度.输出量：向量X的元素是自变量的取值，向量Y的元素是利用（10.15）式、（10.23）式和自适应方法求出初值问题（10.5）在向量X的元素处的数值解，向量H是步长h序列，n是自变量取值的序号，k是自变量的取值.并画出数值解向量Y的图形.function H,X,Y,k,h,P=gaiEuler(funfcn,x0,b,y0,tol)%初始化.pow=1/3; if nargin 5 | isempty(tol), tol = 1.e-6; end;if nargin 6 | isempty(trace), trace = 0; end;x=x0; h=0.0078125*(b-x); y=y0(:);p=128; n=fix(b-x0)/h);H=zeros(p,1); X=zeros(p,1); Y=zeros(p,length(y); k=1; X(k)=x; Y(k,:)=y;% 绘图.clc,x,h,y end% 主循环.while (xx) if x+hbh=b-x; end %计算斜率. fxy=feval(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); endH(k)=h;X(k)=x;Y(k,:)=y; plot(X,Y,mh), gridxlabel(自变量 X), ylabel(因变量 Y)title(用改进的欧拉公式计算dy/dx=f(x,y)，y(x0)=y0在x0,b上的数值解)end%更新步长.if delta=0.0h=min(h*8,0.9*h*(wucha/delta)pow);endendif (x x0=0;y0=1; b=2;tol=1.e-1; Ht,X,Yt,k,h,Pt=odtixing2(funfcn,x0,b,y0,tol)hold onS1 = 8/3*X-29/9+38/9*exp(-3*X),plot(X,S1,b-)hold off hold onH,X,Y,k,h,P=gaiEuler(funfcn,x0,b,y0,tol)hold off legend(用梯形公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解,dy/dx= 8x-3y-7，y(0)=1在0,2上的精确解, 用改进的欧拉公式计算dy/dx=8x-3y-7，y(0)=1在0,2上的数值解),juwY=S1-Y; xiwY=juwY./Y; L=P,S1,juwY,xiwY运行后屏幕显示取精度为，分别用梯形公式和改进的欧拉公式求此解初值问题在区间上的自变量X处数值解Yt, Y和精确解S1及其图形（见图1013），步长H, Y的相对误差xiwY和绝对误差juwY(略).图1013用梯形和改进的欧拉公式求解的图形可以看出改进欧拉公式和梯形公式的结果几乎相同，并且随着的增加累积误差的增长相当可观.10.5 龙格库塔（Runge-Kutta）方法10.5.2 二阶龙格库塔方法的MATLAB程序现提供根据二阶龙格库塔方法的（10.31）式和（10.34）式，编写求解常微分方程初值问题（10.5）的MATLAB程序.用二阶龙格库塔方法求解常微分方程初值问题（10.5）的数值解的MATLAB主程序输入量：funfcn是函数，fun是初值问题的精确解函数， x0和y0是初值b是自变量的最大值，C=c1,c2,a2,b21 是（10.34）式的参数组成的数组, h是步长.输出量：向量X的元素是自变量的取值，向量Y的元素是利用（10.31）式和（10.34）式求出初值问题（10.5）在向量X的元素处的数值解，fxy是在向量X的元素处的精确解，n是自变量取值的序号， wucha(k+1)=norm(Y(k+1)-Y(k)，wch(k)=norm(fxy(k)-Y(k)，并画出数值解和精确解的图形.function k,X,Y,fxy,wch,wucha,P=RK2(funfcn,fun,x0,b,C,y0,h)x=x0; y=y0;p=128; n=fix(b-x0)/h); 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+1 x=x+h;a2=C(3);b21=C(4);c1=C(1);c2=C(2);x1=x+a2*h;k1=feval(funfcn,x,y);y1=y+b21*h*k1;k2=feval(funfcn,x1,y1);fxy(k)=feval(fun,x);y=y+h*(c1*k1+c2*k2);X(k)=x;Y(k,:)=y;k=k+1;plot(X,Y,mh,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(Y(k-1)-Y(k); wch(k)=norm(fxy(k)-Y(k);endX=X(1:k);Y=Y(1:k,:);fxy=fxy(1:k,:);n=1:k; wucha=wucha(1:k,:);wch=wch(1:k,:); P=n,X,Y,fxy,wch,wucha;例10.5.1 用二阶龙格库塔方法求初值问题的数值解，取，并计算与精确解的误差，画出精确解和数值解的图形.解 在MATLAB工作窗口输入下面的程序 x0=0;b=2;C=1/4,3/4,2/3,2/3;y0=0;h=1/4;k,X,Y,fxy,wch,wucha,P=RK2(funfcn,fun,x0,b,C,y0,h)将运行后计算的结果列入（略），画出精确解和数值解的图形（见图1014）.图1014 用二阶龙格库塔方法求初值问题10.5.3 三阶龙格-库塔方法及其MATLAB程序现提供根据三阶龙格-库塔方法（10.38）式和（10.39）式编写求解常微分方程初值问题（10.5）的MATLAB程序.用三阶龙格-库塔方法求解常微分方程初值问题（10.5）的数值解的MATLAB主程序输入量：funfcn是函数，fun是初值问题的精确解函数， x0和y0是初值b是自变量的最大值，C=c1,c2,c3,a2, a3,b21,b31,b32 是（10.38）式的参数组成的数组, h是步长.输出量：向量X的元素是自变量的取值，向量Y的元素是利用（10.38）式和（10.39）式求出初值问题（10.5）在向量X的元素处的数值解，fxy是在向量X的元素处的精确解，n是自变量取值的序号， wucha(k+1)=norm(Y(k+1)-Y(k)，wch(k)= norm(fxy(k)-Y(k)，并画出数值解和精确解的图形.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);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);