第三篇 第十章 常微分方程（组）求解
hold on S1=8/3*x-29/9+38/9*exp(-3*x),plot(X,S1,'b-')
legend('用向前欧拉公式计算dy/dx=8x-3y-7，y(0)=1在[0,2]上的数值
解',' dy/dx =y-x/(3y)，y(0)=1在[0,2]上的精确解')
hold off
juwY=S1-Y; xiwY=juwY./Y;L=[P,S1,juwY,xiwY]
运行后屏幕显示用向前欧拉公式在[0,2]上的自变量X处数值解Y和精确解S1及其图形，步长H，Y的相对误差xiwY和绝对误差juwY。
10.3.3 向后欧拉公式及其误差分析
y(xn?1)?y(xn) 如果在用差商代替(10.5)方程的导数y?后，f(x,y)中的xh取小区间[xn,xn?1]的右端点xn?1，再同样地以yn,yn?1分别代替y(xn)，y(xn?1)，就得到
yn?1?yn?hf(xn?1,yn?1),n?0,1,2,... (10.15)
称为向后欧拉公式，它也可看作(10.9)式中积分用右端点矩形公式近似的结果。但是(10.15)式右端的yn?1未知，所以称隐式公式，无法用它直接计算yn?1。
隐式公式(10.15)通常用迭代发计算，即先有向前欧拉公式产生初值
0yn,2,... (10.16) ?1?yn?hf(xn,yn),n?0,1再按下式进行迭代
(k?1)(k)yn,2,... (10.17) ?1?yn?hf(xn?1,yn?1),n?0,1若序列yn?1??(k)??k?0收敛，则
(k)ynlim?1?yn?1. k??? 为了估计 (10.15)式的局部误差。将右端的f(x,y(x))在精确点(xn?1,y(xn?1))展开，并注意到yn?y(xn)，有
yn?1?y(xn)?h?f(xn?1,y(xn?1))?fy(xn?1,?n?1)?yn?1?y(xn?1)?? (10.18)
其中?n?1在y(xn?1)与yn?1之间，再将
f(xn?1,y(xn?1))?y?(xn?1)?y?(xn)?hy??(xn)?O(h2) (10.19)
代入 (10.18)后，与 (10.12)式相减，得
1y(xn?1)?yn?1??h2y??(xn)?O(h3)?hfy(xn?1,?n?1)?y(xn?1)?yn?1?, (10.20)
2于是
1?h2y??(xn)?O(h3)1y(xn?1)?yn?1?2??h2y??(xn)?O(h3)?O(h2). (10.21)
1?hfy(xn?1,?n?1)2与向前欧拉公式的误差估计(10.13)式相比，我们注意到它们的h阶项数值相等，而符号
相反，但二者都只有1阶精度。
210.3.4 向后欧拉方法的MATLAB程序
用向后欧拉方法求解常微分方程初值问题（10.5）的数值解的MATLAB主程序
输入量：funfcn是函数f(x,y)，x0和y0是初值问题y(x0)?y0,b是自变量x的最大值，h是步长，tol是精度。
输出量：向量X的元素是由自变量x的取值组成，数组Y的元素是利用(10.16)式和(10.17)式求出常微分方程初值问题(10.5)在X的元素处的数值解，n是自变量x取值的个数。并画出数值解向量Y的图形。
function [X,Y,n,P]=Heuler1(funfcn,x0,b,y0,h,tol)
185．
第三篇 第十章 常微分方程（组）求解
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*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
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]上的数
值解')
end
X=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1; P=[n',X,Y]
例10.3.6 用向后欧拉公式求解区间[0,2]上的初值问题
dy??3y?8x?7, dxy(0)?1的数值解，取步长h?0.05，并与精确解作比较，在同一个坐标系中作出图形.然后再取h?0.01，观察数值解与精确解误差的变化，说明h与误差的关系.
解 (1)建立并保存名为funfcn.m的M文件函数。 function f=function(x,y) f=8*x-3*y-7;
(2)建立并保存名为Heuler1.m的M文件函数。 (3)输入程序
>>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 on
S2= 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 off
juwY1=S2-Y1;
xiwY1=juwY1./Y1;
L=[P1,S2,juwY1,xiwY1] subplot(2,1,2) h=0.05;
[X,Y,n,P]=Heuler1(@funfcn,x0,b,y0,h,tol) hold on
S1 = 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 off
186．
第三篇 第十章 常微分方程（组）求解
juwY=S1-Y; xiwY=juwY./Y;
L=[P,S1,juwY,xiwY]
运行后屏幕显示用向后欧拉公式计算此初值问题在[0,2]上的自变量X处数值解Y和精确解S1及其图形，步长H, Y的相对误差xiwY和绝对误差juwY(略) . 习题三：
?dy??1?y?x,0?x?1,1，分别用向前欧拉方法和向后欧拉方法求初值问题?dx的数值解，
??y(0)?1分别取h?0.075,0.0075，并计算误差，画出精确解和数值解的图形。
?dy22??1?y?x,0?x?1,2，用向前欧拉公式（10.8）求解初值问题?dx，取n?10,50，
??y(0)?0并将计算结果与精确解作比较，写出在每个子区间[xk,xk?1]上的局部截断误差公式，画
出精确解和数值解在区间[0，1]上的图形。 3，用向后欧拉公式求y(x)??x0e?2tdt在x?0.5,1.0,1.5,2.0处的近似值。
24，分别用向前欧拉公式和一阶泰勒多项式逼近法及其函数double求解初值问题
?dy??1?ycosx,0?x?1,，取n?40,80时，并将计算结果与精确解作比较，写出在每?dx??y(0)?0个子区间[xk,xk?1]上的局部截断误差公式，画出精确解和数值解在区间[0，1]上的图形，
并说明用函数double计算的结果是否有实用价值。
5，用自适应向前欧拉方法求解在区间[0,2]上的初值问题
dy??y?x2?x,y(0)?0，取精度10?1并与精确解作比较，并在同一个坐标系中作出图dx形。
10.4 改进的欧拉方法的MATLAB程序
10.4.1 梯形公式及其误差估计
将向前欧拉公式(10.10)和向后欧拉公式和(10.15)加以平均，得到
hyn?1?yn??f(xn,yn)?f(xn?1,yn?1)?,n?0,1,2,... (10.22)
2它相当于(10.9)式：y(xn?1)?yn?1??xn?1xnf?x,y(x)?dx(n?0,1,2,...)右端的积分用
求面积的梯形公式近似，再以yn，yn?1分别代替y(xn)，y(xn?1)的结果，(10.22)式称为梯形公式。它也要用迭代方法求解，即由(10.16)式
0yn,2,...得到初值后按照下失进行迭代 ?1?yn?hf(xn,yn),n?0,1h(k?1)yn?y??f(xn,yn)?f(xn?1,yn?1)?,k?0,1,2,...,n?0,1,2,... (10.2 3)?1n2因为式(10.13)
123 y(xn?1)?yn?1?hy??(xn)?O(h)
2和(10.21)式
187．
第三篇 第十章 常微分方程（组）求解
y(xn?1)?yn?1??212hy??(xn)?O(h3) 2的右端h阶象正好抵消，所以梯形公式的局部截断误差为O(h3)，具有2阶精度。
10.4.2 梯形公式的两种MATLAB程序
(一) 梯形公式的MATLAB程序
用梯形公式求解常微分方程初值问题（10.5）的数值解的MATLAB主程序1
输入量：funfcn是函数f(x,y)，x0和y0是初值问题y(x0)?y0,b是自变量x的最大值，h是步长，tol是精度。
输出量：向量X的元素是由自变量x的取值组成，数组Y的元素是利用(10.16)式、(10.23)式和自适应方法求出常微分方程初值问题(10.5)在X的元素处的数值解，向量H是步长h序列，n是自变量x取值的序号。并画出数值解向量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 Wu>tol p=Y1(i,:),
fxip=feval(funfcn,X(i),p);
Y1(i,:)=y0+h*(fx0y0+fxip)/2,P1=Y1(i,:), Y(i,:)=p1; 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]上的数值
解')
end
X=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1; P=[n',X,Y]
例10.4.1 用梯形公式求解区间[0,2]上的初值问题
?1dy??3y?8x?7,y(0)?1，dx取步长h?0.05，精度为10，并与精确解作比较，在同一个坐标系中画出图形.
解 (1)建立并保存名为funfcn.m的M文件函数。 function f=function(x,y) f=8*x-3*y-7;
(2)建立并保存名为odtixing1.m的M文件函数。
188．
第三篇 第十章 常微分方程（组）求解
(3)输入程序
>>S1=dsolve('Dy=8*x-3*y-7','y(0)=1','x') >> x0=0;y0=1; b=2; tol=0.1; h=0.05;
[X,Yt,n,Pt]=odtixing1(@funfcn,x0,b,y0,h,tol) hold on
S1=8/3*X-29/9+38/9*exp(-3*X); plot(X,S1,'b-'), hold off
legend('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] 运行后屏幕显示取精度为10，分别用梯形公式和向前欧拉公式求解此初值问题在区间[0,2]上的自变量X处数值解Yi(i=t,q)和精确解S1，步长H, Yi的相对误差xiwYi和绝对误差juwYi (略) 及其数值解和精确解的图形.
(二) 自适应梯形公式的MATLAB程序
用自适应梯形公式求解常微分方程初值问题（10.5）的数值解的MATLAB主程序 输入量：funfcn是函数f(x,y)，x0和y0是初值问题y(x0)?y0,b是自变量x的最
大值，tol是精度。
输出量：向量X的元素是由自变量x的取值组成，数组Y的元素是利用(10.16)式、(10.23)式和自适应方法求出常微分方程初值问题(10.5)在X的元素处的数值解，向量H是步长h序列，n是自变量x取值的序号。并画出数值解向量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 (x
%计算斜率.
fxy=feval(funfcn,x,y); fxy=fxy(:); %计算误差，设定可接受误差. delta=norm(h*fxy,'inf');
wucha=tol*max(norm(y,'inf'),1.0);
% 当误差可接受时重写解. if delta<=wucha
x=x+h; y1=y+h*fxy; fxy1=feval(funfcn,x,y1); fxy=fxy(:); y2=y+h*fxy1; y=(y1+y2)/2; k=k+1; if k>length(X)
X=[X;zeros(p,1)]; Y=[Y;zeros(p,length(y))]; H=[H;zeros(p,1)]; end
189．
?1
百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库第十章 常微分方程(组)求解(4)在线全文阅读。
相关推荐：