第10章 - Matlab与数值分析

这里n是等分段数，h?(b?a)/n是步长，节点为xi?a?kh,k?0,1,?,n．

2例10.27 用逐次分半的复化梯形公式计算I??4xsinxdx10,精确到. ?0% 逐次分半的复化梯形积分公式

fun=inline('x.*sin(x))');% 输入函数 a=0;b=2;% 输入区间

h=b-a; % 步长预设为区间长 e=1e-4; % 要求精度 % 以下计算积分值

T=(feval(fun,a)+feval(fun,b))/2*h n=1;

r=1; % 当前误差预设为1

while r>e % 当前误差小于要求精度时退出计算 h=h/2;n=2*n; s=0;

for k=2:2:n

s=s+feval(fun,a+(k-1)*h); end T0=T;

T=T0/2+s*h % 逐次分半公式 r=abs(T0-T); % 计算当前误差 end T

运行结果： T = 1.7416

3． 复化Simpson公式

mm?1h??I??f(a)?f(b)?2?f(x2k)?4?f(x2k?1)?

3?k?0k?1?这里２m是等分段数，h?b?a是步长，节点为xi?a?kh,k?0,1,?,2m． 2m2例10.28 用分10段的复化Simpson公式计算I?%quad_4 复化Simpson积分公式 fun=inline('x.*sin(x)'); % 输入函数 a=0;b=2;% 输入区间 m=10;% 输入等分段数 h=(b-a)/m/2;% 计算步长 % 以下计算积分值

S=feval(fun,a)+feval(fun,b);

218

?4xsinxdx10,精确到. ?0for k=2:m

S=S+2*feval(fun,a+2*(k-1)*h); end

for k=1:m

S=S+4*feval(fun,a+(2*k-1)*h); end S=S*h/3 运行结果： S = 1.7416

10.8.3 Romberg积分计算

以下编制的是四列形式的Romberg公式，

其中第一个序列是T序列，它由逐次分半的梯形公式形成，收敛较慢； 第二个序列是S序列，由T序列的线性组合Si?1?4Ti?1?Ti形成； 316Si?1?Si形成；

1564Ci?1?Ci形成．

63第三个序列是C序列，由S序列的线性组合Ci?1?第四个序列是R序列，由C序列的线性组合Ri?1?2例10.29 用Romberg公式计算I??xsinxdx，当Ri?1?Ri0???10?5时退出．

编程计算：

% Romberg积分公式

fun=inline('x.*sin(x)'); % 输入函数 a=0;b=2; % 输入区间 e=1e-6; % 输入精度要求 % 以下是Romberg过程 h=b-a;n=1;

r=1; % 当前误差预设为1 T=[];S=T;C=T;R=T;

T(1)=(feval(fun,a)+feval(fun,b))/2*h; i=2;

while r>e % 当前误差小于要求精度时退出计算 h=h/2;n=2*n; s=0;

for k=2:2:n

s=s+feval(fun,a+(k-1)*h); end

T(i)=T(i-1)/2+s*h; % 形成T序列

219

S(i)=(4*T(i)-T(i-1))/3; % 形成S序列 if i>=3

C(i)=(16*S(i)-S(i-1))/15; % 形成C序列 end if i>=4

R(i)=(64*C(i)-C(i-1))/63; % 形成R序列 end if i>4

r=abs(R(i)-R(i-1)); % 重新计算当前误差 end i=i+1; end

A=[T',S',C',R'] % 得到结果 运行结果： A =

1.8186 0 0 0 1.7508 1.7282 0 0 1.7434 1.7409 1.7417 0 1.7420 1.7415 1.7416 1.7416

1.7417 1.7416 1.7416 1.7416

10.8.4 Matlab中的积分公式

在Matlab中有一批积分命令，直接使用非常简单方便． 1. 复化梯形公式：Trapz( )

格式：I=trapz(x,y) 其中x是离散的数据点，y是对应的函数值． 同上例用n=10的复化梯形公式计算积分值： >> fun='x.*sin(x)'; >> x=0:1/10:2; >> y=eval(fun)； >> I=trapz(x,y) I = 1.7417

2. 复化Simpson公式：quad( ) 格式：I=quadl(fun,a,b,e) e是积分的精度要求，缺省值是10-6．它采用的是自适应的Simpson公式，即对Simpson公式逐次对分，直到满足精度要求为止． 同上例用自适应的复化Simpson公式计算积分值，精度要求10-6： >> fun='x.*sin(x)'; >> I=quad(fun,0,2) I = 1.7416

3. 复化Lobatto公式：quadl( )

格式：I=quadl(fun,a,b,e) 它采用的是自适应的复化Lobatto公式，Lobatto公式是一类闭式Gauss型积分公式，代数精度较高．用法同quad( )． >> I=quadl(fun,0,2) I = 1.7416.

220

4. 复化的8阶Newton－Cotes公式：quad8( )

格式：I=quad8(fun,a,b,e) 它采用的是自适应的复化8阶Newton－Cotes公式，代数精度较高．用法同quad( )． >> I=quad8(fun,0,2) I = 1.7416

5. 二重积分和三重积分公式 二重积分公式：dblquad( ); 三重积分公式：traplequad( );

这两个公式参数较多，使用时请参阅有关资料．仅以dblquad( )举例说明： 例：计算I???Dy?x2dxdy, D??(x,y)|?1?x?1,0?y?2?，精度要求10-6

>> fun='sqrt(abs(y-x.^2))'； >> I=dblquad(fun,-1,1,0,2) I = 3.2375

§10.9 常微分方程初值问题数值解

?y??f(x,y) a?x?b ?y(a)?y0?(x) a?求 y?yx? b称为常微分方程初值问题，有多种数值方法求解．

10.9.1 单步法

１．Euler方法

?y0?y(a) ?y?y?hf(x,y) k?0,1,2,?kkk?k?1例10.30 取h=0.1, 用Euler方法计算:

2x?y?y? 0?x?1 ?y ??y(0)?1?%Euler方法

%输入函数f(x,y)和求解区间[a,b] fun='y-2*x./y'; a=0;b=1; h=0.1;

n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1); %Euler方法 X(1)=a;Y(1)=1; for i=1:n

x=X(i);y=Y(i);

221