哈工大-数值分析上机实验报告

时间：2020-09-01 08:12:54　来源：勤学考试网本文已影响人

实验报告

题目：非线性程求解

摘要：非线性程的解析解通常很难给出，因此线性程的数值解法就尤为重要。本实验采用两种常见的求解法二分法和 Newton法及改进的Newton法。

前言：(目的和意义)掌握二分法与 Newton法的基本原理和应用。

　数学原理：

对丁一个非线性程的数值解法很多。在此介绍两种最常见的法：二分法和Newton法'

对丁二分法，其数学实质就是说对丁给定的待求解的程 f(x),其在［a,b］上连续，

f(a)f(b)<0，且f(x)在［a,b］仅有一个实根x,取区间中点c,若，贝U c恰为其根，否则根据 f(a)f(c)<0是否成立判断根在区间［a,c］和［c,b］中的哪一个，从而得出新区问，仍称为［a,b］。

　重复运行计算，直至满足精度为止。这就是二分法的计算思想。

Newton法通常预先要给出一个猜测初值 x°,然后根据其迭代公式

Xk 1 xk

f (xk)

产生逼近解x的迭代数列{xk},这就是Newton法的思想。当x°接近x时收敛很快，但是当xo选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为

f (xk)

xk 1 xk r

f (xk)

其中r为要求的程的根的重数，这就是改进的 Newton法，当求解已知重数的程的根时,

在同种条件下其收敛速度要比 Newton法快的多。

程序设计：

本实验采用Matlab的M文件编写。其中待求解的程写成 function的式，如下

function y=f(x);

y=-x*x-sin(x);

写成如上形式即可，下面给出主程序。

二分法源程序：

clear

%%%定求解区问

b=1.5;

a=0;

%%%差

R=1;

k=0;以迭代次数初值

while (R>5e-6);

c=(a+b)/2;

if f12(a)*f12(c)>0;

a=c;

else

b=c;

end

R=b-a;%求出误差

k=k+1;

end

x=c%合出解

Newton法及改进的 Newton 法源程序：

clear

%%%%俞入函数

f=input('请输入需要求解函数>>','s')

%%%解f(x)的导数

df=diff(f);

%%%进常数或重根数

miu=2;

%%%始值x0

x0=input('input initial value x0>>');

k=0;%!代次数

max=100 ;%!大迭代次数

R=eval(subs(f,'x0','x'));咻解f(x0)，以确定初值x0时否就是解

while (abs(R)>1e-8)

x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x'));

R=x1-x0;

x0=x1;

k=k+1;

if (eval(subs(f,'x0','x'))<1e-10);

break

end

if k>max;%fc果迭代次数大丁给定值，认为迭代不收敛，重新输入初值

ss=input('maybe result is error,choose a new x0,y/n?>>','s');

if strcmp(ss,'y')

x0=input('input initial value x0>>');

k=0;

else

break

end

k;以给出迭代次数

x=x0；以%合出解

结果分析和讨论：

用二分法计算程sin x ： 0在［1 , 2］的根。(5* 10 6，下同) 计算结果为

x= 1.523 ;

f(x)= -3.4311e-007 ;

k=18 ;

由f(x)知结果满足要求，但迭代次数比较多，法收敛速度比较慢。

用二分法计算程x3 x 1 0在［1, 1.5］的根。

计算结果为

x= 1.180 ;

f(x)= 2.4815e-006 ;

k=17 ;

由f(x)知结果满足要求，但迭代次数还是比较多。

用Newton法求解下歹U程

xe 1 0 x0=0.5;

计算结果为

x= 0.978 ;

f(x)= 2.0313e-016 ;

k=4 ;

由f(x)知结果满足要求，而且乂迭代次数只有 4次看出收敛速度很快。

x3 x 1 0 x°=1 ;

(x 1)2 (2x 1) 0 x°=0.45, x 0=0.65 ;

当x0=0.45时，计算结果为

x= 0.983 ;

f(x)= -8.4584e-014 ;

k=4; 由f(x)知结果满足要求，而且乂迭代次数只有 4次看出收敛速度很快，实际上该程确实有

真解x=0.5 o

当xo=0.65时，计算结果为 x= 0.000 ; f(x)=0; k=9 ; 由f(x)知结果满足要求，实际上该程确实有真解 x=0.5,但迭代次数增多，实际上当取 x0>

0.68时，xE ,就变成了程的另一个解，这说明 Newton法收敛与初值很有关系，有的时候甚至可能不收敛。

用改进的Newton法求解，有2重根，取 2

(x 1)2(2x 1) 0 河=0.55;并与3.中的c)比较结果。

当x0=0.55时，程序死循环，无法计算，也就是说不收敛。改 1.5时，结果收敛为 x=0.286 ; f(x)=4.1127e-007 ; k=16 ; 显然这个结果不是很好，而且也不是收敛至程的 2重根上。

当xo=0.85时，结果收敛为 x= 1.489 ; f(x)= 2.8737e-; k=4 ; 这次达到了预期的结果，这说明初值的选取很重要，直接关系到法的收敛性，实际上直接

用Newton法，在给定同样的条件和精度要求下，可得其迭代次数 k=15 ,这说明改进后

的Newton法法速度确实比较快。

　结论：对丁二分法，只要能够保证在给定的区间有根，使能够收敛的，当时收敛的速度和给

定的区间有关，二且总体上来说速度比较慢。 Newton法，收敛速度要比二分法快，但是

最终其收敛的结果与初值的选取有关，初值不同，收敛的结果也可能不一样，也就是结果

可能不时预期需要得结果。改进的 Newton法求解重根问题时，如果初值不当，可能会不收敛，这一点非常重要，当然初值合适，相同情况下其速度要比 Newton法快得多。

实验报告二

题目： Gauss歹U主元消去法

摘要：求解线性程组的法很多，主要分为直接法和间接法。本实验运用直接法的 Guass

消去法，并采用选主元的法对程组进行求解。

前言：(目的和意义)

学习Gauss消去法的原理。

了解列主元的意义。

确定什么时候系数阵要选主元数学原理：

由丁一般线性程在使用 Gauss消去法求解时，从求解的过程中可以看到，若a$ 1)=0 ,

则必须进行行交换，才能使消去过程进行下去。有的时候即使 a；：1〉 0,但是其绝对值非

常小，由丁机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。因此有必要进行列主元技术，以最大可能的消除这种现象。这一技术要寻找行 r,使得

闻 1)| max a" 1)

I k

并将第r行和第k行的元素进行交换，以使得当前的a：： 1)的数值比0要大的多。这种列主元的消去法的主要步骤如下：

消元过程

对k=1,2, ；n-1,进行如下步骤。

选主元，记

I ark I max 孤

i k

若|ark I很小，这说明程的系数矩阵重病态，给出警告，提示结果可能不对。

交换增广阵A的r, k两行的元素。

a「j &灯 (j=k, ?-；n+1)

计算消元

a0 a0 aikakj/akk (i=k+1,? ?,n; j=k+1, …,n+1)

回代过程

对k= n, n-1, ?-,1,进行如下计算

xk (ak,n 1 akj xj / akk )

至此，完成了整个程组的求解。

程序设计：

本实验采用Matlab的M文件编写。

Gauss消去法源程序：

clear

a=input('输入系数阵：>>\n')

b=input('输入歹U阵 b: >>\n')

n=length(b);

A=[a b]

x=zeros(n,1);

%% %数主体

for k=1:n-1;

%%%否进行主元选取

if abs(A(k,k))<yipusilong; %事先给定的认为有必要选主元的小数 yzhuyuan=1;

else yzhuyuan=0;

end

if yzhuyuan;

%%%^ 元

t=A(k,k);

for r=k+1:n;

if abs(A(r,k))>abs(t)

p=r;

else p=k;

end

%%%换元素

if p~=k;

for q=k:n+1;

s=A(k,q);

A(k,q)=A(p,q);

A(p,q)=s;

end

%%%断系数矩阵是否奇异或病态非常重

if abs(A(k,k))< yipusilong

disp(矩阵奇异，解可能不正确’)

end

%%%馈消元，得三角阵

for r=k+1:n;

m=A(r,k)/A(k,k);

for q=k:n+1;

A(r,q)=A(r,q)-A(k,q)*m;

end

%%%%牟 x

x(n)=A(n,n+1)/A(n,n);

for k=n-1:-1:1;

s=0;

for r=k+1:n;

s=s+A(k,r)*x(r);

end

t=(A(k,n+1)-s)

x(k)=(A(k,n+1)-s)/A(k,k)

end

结果分析和讨论：

x 22

例：求解程

y 34。其中为一小数，当

10 5,10 10,10 14,10 20 时,分

z 10

别米用列主元和不列主元的

Gauss消去法求解，并比较分果。

记Emax为求出的解代入程后的最大误差，按要求，计算结果如下：

当 10 5时，不选主元和选主元的计算结果如下，其中前一列为不选主元结果，后

一歹0为选主元结果，下同。

0.391

0.651

2.972

2.163

2.451

2.721

Emax= 9.2624e- , 0

此时，由丁不是很小，机器误差就不是很大，由 Emax可以看出不选主元的计算结果

精度还可以，因此此时可以考虑不选主元以减少计算量。

当 10 -15，所以此时的 10 20就认为是0,故出现了错误现象。而选主元时则没有这种现象，

-15，所以此时的 10 20就认为是0,故出现了错误现象。而选主元时则没有这种现象，

而且由Emax可以看出选主元时的结果应该是精确解。

结论：

1.877

0.348

1.807

2.174

3.731

2.609

Emax= 2.4668e-005 , 0

此时由Emax可以看出不选主元的计算精度就不好了，误差开始增大。

当 10 米用Gauss消去法时，如果在消兀时对角线上的兀素始终较大（假如大丁 10 ）,那

米用Gauss消去法时，如果在消兀时对角线上的兀素始终较大（假如大丁 10 ）,那

么本法不需要进行列主元计算，计算结果一般就可以达到要求，否则必须进行列主元这一

步，以减少机器误差带来的影响，使法得出的结果正确。

1.020

1.000

1.666

2.000

3.111

0000

Emax= 0.503,

此时由Emax可以看出，不选主元的结果应该可以说是不正确了，这是由机器误差引起

当 10 20时，不选主元和选主元的计算结果如下

NaN

Emax=NaN, 0

不选主元时，程序报错：

Warning: Divide by zero.。这是因为机器计算的最小精度为

的。

实验报告三

题目：Rung现象产生和克服

摘要：由丁高次多项式插值不收敛，会产生 Runge现象，本实验在给出具体的实例后，

采用分段线性插值和三次样条插值的法有效的克服了这一现象，而且还取的很好的插值效

果。

前言：(目的和意义)

深刻认识多项式插值的缺点。

明确插值的不收敛性怎样克服。

明确精度与节点和插值法的关系。

数学原理：

在给定n+1个节点和相应的函数值以后构造 n次的Lagrange插值多项式，实验结果表明(见后面的图)这种多项式并不是随着次数的升高对函数的逼近越来越好，这种现象

就是Rung现象。

解决Rung现象的法通常有分段线性插值、三次样条插值等法。

分段线性插值：

设在区间［a, b］上，给定n+1个插值节点

a=x0<xi< <xn=b

和相应的函数值y。，yi，…，yn,,求作一个插值函数 (x),具有如下性质：

(x。yj, j=0, 1,…，n。

(x)在每个区间［xi, xj］上是线性连续函数。则插值函数 (x)称为区间［a, b］上对应n

个数据点的分段线性插值函数。

三次样条插值：

给定区问［a, b］一个分划

/: a=xo<x1< <xN=b

若函数S(x)^足下列条件：

S(x)在每个区间［xi, x］上是不高丁 3次的多项式。

S(x)及其2阶导数在［a, b］上连续。则称S(x)ffi关丁分划/的三次样条函数。

　程序设计：

本实验采用Matlab的M文件编写。其中待插值的程写成 function的式，如下

function y=f(x); y=1/ (1+25*x*x);

写成如上形式即可，下面给出主程序

Lagrange插值源程序：

n=input('将区间分为的等份数输入：\n');

s=[-1+2/n*[0:n]]; %%%定的定点，Rf为给定的函数 x=-1:0.01:1;

f=0;

for q=1:n+1;

l=1;%求插值基函数

for k=1:n+1;

if k~=q;

l=l*(x-s(k))./(s(q)-s(k));

else

l=l;

end

f=f+Rf(s(q))*l; %求插值函数

end

plot(x,f,'r') %乍出插值函数曲线

grid on

hold on

分段线性插值源程序

clear

n=input('将区间分为的等份数输入：\n');

s=[-1+2/n*[0:n]]; %%%定的定点，Rf为给定的函数 m=0;

hh=0.001;

for x=-1:hh:1;

ff=0;

for k=1:n+1; %%%插值基函数

switch k

case 1

if x<=s (2);

l=(x-s(2))./(s(1)-s(2));

else

l=0;

end

case n+1

if x>s(n);

l=(x-s(n))./(s(n+1)-s(n));

else

l=0;

end

otherwise

if x>=s(k-1)&x<=s(k);

l=(x-s(k-1))./(s(k)-s(k-1));

else if x>=s(k)&x<=s(k+1);

l=(x-s(k+1))./(s(k)-s(k+1));

else

l=0;

end

ff=ff+Rf(s(k))*l; %%:插值函数值

end

m=m+1;

f(m)=ff;

end

%%%出曲线

x=-1:hh:1;

plot(x,f,'r');

grid on

hold on

三次样条插值源程序：(采用第一边界条件)

clear

n=input('将区间分为的等份数输入：\n');

%遍值区间

a=-1;

b=1;

hh=0.001 ;%U图的步长

s=[a+(b-a)/n*[0:n]]; %%%定的定点，Rf为给定的函数

%%%W边界条件 Rf"(-1),Rf"(1)

vn5000*qd+25*a*a)>3&0d+25*a*a)>4 forknm

h(kHs(k+vs(k)」

end

fork"—km——k;粤任洲密回*-amudamiu _a(kHh(k-M)、(h(k-M)+h(k))」 miu(kHq'a(k)」

end

汶汶盅鬟次密盗番A

for kn_km_k-

for pn_k5—k-

swioh p

case k

akphn

case k—_k

A(k-pHmiu(p+q)； case k+_k

A(k-P)*(pa;

otherwise

akphp

end

汶汶d番

for kn_km_k-

swioh k

case—k

d(kH6*f2c(『s(k) s(k+q) s(k+2)vmiu(k)*< case rk_k

d(kH6*f2c(『s(k) s(k+q) s(k+2)D'a(k)*< otherwise

d(kH6*f2c(『s(k) s(k+q) S(k+2)M end

end

汶汶卷海M番

M￡d」

4怜漆驾

M=[v;M;v];

%%%%

m=0; f=0; for x=a:hh:b; if x==a;

p=1;

else

p=ceil((x-s(1))/((b-a)/n));

end ff1=0; ff2=0; ff3=0; ff4=0; m=m+1; ff1=1/h(p)*(s(p+1)-x)A3*M(p)/6; ff2=1/h(p)*(x-s(p))A3*M(p+1)/6; ff3=((Rf(s(p+1))-Rf(s(p)))/h(p)-h(p)*(M(p+1)-M(p))/6)*(x-s(p)); ff4=Rf(s(p))-M(p)*h(p)*h(p)/6; f(m)=ff1+ff2+ff3+ff4 ; end %%%出插值图形 x=a:hh:b; plot(x,f,'k') hold on grid on

结果分析和讨论：

本实验米用函数f (x) 2进行数值插值，插值区间为[-1, 1],给定节点为

1 25x2

xj=-1 +jh,h=0.1, j=0,…,n。下面分别给出Lagrange插值，三次样条插值，线性插值的函数曲线和数据表。图中只标出 Lagrange插值的十次多项式的曲线，其它曲线没有标出，从数据表中可以看出具体的误差。

表中，Lw（x）为Lagrang e插值的10次多项式，Si0（x）, S40（x）分别代表n=10 , 40的三次样条插值函数，Xi0（x）, X40（x）分别代表n=10 , 40的线性分段插值函数。

f（x）

Li0（x）

S10（x）

-1.000

0.154

-0.000

0.239

1.920

0.040

0.239

0.910

0.239

-0.000

0.941

1.926

0.458

0.941

0.665

0.941

-0.000

0.344

0.982

0.979

0.344

0.421

0.344

-0.000

0.176

-0.000

0.378

-0.674

0.744

0.378

0.882

0.378

-0.000

0.321

-0.250

0.866

0.321

0.588

0.321

-0.000

0.649

-0.418

0.849

0.649

0.294

0.649

-0.000

0.000

-0.000

0.788

0.257

0.713

0.788

0.000

0.788

-0.000

0.276

0.103

0.730

0.276

0.000

0.276

-0.000

0.825

0.267

0.883

0.825

0.000

0.825

-0.000

0.000

-0.000

0.385

0.376

0.464

0.385

0.000

0.385

-0.000

0.231

0.080

0.860

0.231

0.000

0.231

-0.000

0.902

0.789

0.327

0.902

0.000

0.902

-0.000

0.000

-0.000

0.000

0.340

0.431

0.000

-0.000

0.000

0.890

0.828

0.000

-0.000

0.824

0.073

0.810

0.824

0.000

0.824

0 1.000 1.000 1.000 1.000 1.000

1.000

0.000

0.824

0.073

0.810

0.824

0.000

0.824

0.000

0.890

0.828

0.000

0.340

0.431

0.000

0.902

0.789

0.327

0.902

0.000

0.902

0.000

0.231

0.080

0.860

0.231

0.000

0.231

0.000

0.385

0.376

0.464

0.385

0.000

0.385

0.000

0.825

0.267

0.883

0.825

0.000

0.825

0.000

0.276

0.103

0.730

0.276

0.000

0.276

S40（x） Xi0（x） X40（x）

0.000

0.788

0.257

0.713

0.788

0.000

0.788

0.000

0.649

-0.418

0.849

0.649

0.294

0.649

0.000

0.321

-0.250

0.866

0.321

0.588

0.321

0.000

0.378

-0.674

0.744

0.378

0.882

0.378

0.000

0.176

0.000

0.344

0.982

0.979

0.344

0.421

0.344

0.000

0.941

1.926

0.458

0.941

0.665

0.941

0.000

0.239

1.920

0.040

0.239

0.910

0.239

1.000

0.154

从以上结果可以看到，用三次样条插值和线性分段插值，不会出现多项式插值是出现

的Runge现象，插值效果明显提高。进一步说，为了提高插值精度，用三次样条插值和线性分段插值是可以增加插值节点的办法来满足要求，而用多项式插值函数时，节点数的

增加必然会使多项式的次数增加，这样会引起数值不稳定，所以说这两种插值要比多项式插值好的多。而且在给定节点数的条件下，三次样条插值的精度要优丁线性分段插值，曲

线的光滑性也要好一些。

实验报告四

题目：多项式最小二乘法摘要：对丁具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观

察和测量给出离散的一些点，再来求出具体的函数解析式。

　乂因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。最小二乘法是求解这一问题的很好的

法，本实验运用这一法实现对给定数据的拟合。

前言：（目的和意义）

学习使用最小二成法的原理

了解法程的特性数学原理:

对丁给定的测量数据（x,fi）（i=1,2，…，n）,设函数分布为

y(x)m

y(x)

aj j(x)

j 0

特别的，取j（x）为多项式

j(x) x

j(x) xj

(j=0, 1,…，m)

则根据最小二乘法原理，可以构造泛函

H (a°,ai,,

H (a°,ai,

,am )

(fi

i 1

aj j(xi))

j 0

(k=0, 1,…，m)

则可以得到法程

0)1)0)1)0)1)a。

a。

a〔

(f,

（

求该解程组，则可以得到解

a0,a1,

（f,

,am,因此可得到数据的最小二乘解

f(x) aj j(x)

程序设计：

本实验采用Matlab的M文件编写。其中多项式函数 j xj写成function的式，如

下

function y=fai(x,j)

y=i;

for i=1:j

y=x.*y;

end

写成如上形式即可，下面给出主程序。

多项式最小二乘法源程序

clear

%%%定测量数据点(s,f)

s=[3 4 5 6 7 8 9];

f=[2.01 2.98 3.50 5.02 5.47 6.02 7.05];

%%%算给定的数据点的数目

n=length(f);

%%%定需要拟合的数据的最高次多项式的次数

m=10;

%%%序主体

for k=0:m;

g=zeros(1,m+1);

for j=0:m;

t=0;

for i=1:n;咐算积(fai(si),fai(si))

t=t+fai(s(i),j)*fai(s(i),k);

end

g(j+1)=t;

end

A(k+1,:)=g;姒程的系数矩阵

t=0;

for i=1:n;咐算积(f(si),fai(si))

t=t+f(i)*fai(s(i),k);

end

b(k+1,1)=t;

end

a=A\b咻出多项式系数 x=[s(1):0.01:s(n)]';

y=0;

for i=0:m;

y=y+a(i+1)*fai(x,i);

end

plot(x,y) %乍出拟合成的多项式的曲线

grid on

hold on plot(s,f,'rx') %E上图中标记给定的点

结果分析和讨论：

例用最小二乘法处理下面的实验数据

2.01

2.98

3.50

5.02

5.47

6.02

7.05

并作出f (x)的近似分布图。

分别采用一次，二次和五次多项式来拟合数据得到相应的拟合多项式为:

y1=-0.38643+0.82750x ;

y2=-1.03024+1.06893x-0.02012x 2;

y5=-50.75309+51.53527x-19.65947x 2+3.66585x 3-0.32886x 4+0.01137x5;

分别作出它们的曲线图，图中点划线为 y1曲线，实线为y2曲线，虚线为y5曲线。’x'为给定的数据点。从图中可以看出并不是多项式次数越高越好，次数高了，曲线越能给定点

实验报告五

题目： Romberg积分法

摘要：对丁实际的工程积分问题，彳艮难应用 Newton-Leibnitz公式去求解。因此应用数值法进行求解积分问题已经有着很广泛的应用，本文基丁 Romberg积分法来解决一类积分

|可题。

前言：(目的和意义)

理解和掌握Romberg积分法的原理；

学会使用Romberg积分法；

明确Romberg积分法的收敛速度及应用时容易出现的问题。

数学原理:

考虑积分I(f) f (x)dx，欲求其近似值，通吊有复化的梯形公式、 Simpsion公式

和Cotes公式。但是给定一个精度，这些公式达到要求的速度很缓慢。如提高收敛速度，

自然是人们极为关心的课题。为此，记「k为将区间［a,b］进行2k等分的复化的梯形公式计

算结果，记T2,k为将区间［a,b］进行2k等分的复化的Simpsion公式计算结果，记 T3,k为将区问［a,b］进行2k等分的复化的Cotes公式计算结果。根据 Richardson外推加速法，可以得到收敛速度较快的Romberg积分法。其具体的计算公式为：

准备初值，计算

a b

Ti,i 亍[f(a) f(b)]

2.按梯形公式的递推关系，计算

2k 1 1

"孔蜡 f(a -r^(i 0.5))

2 2 i 0 2

3.按Romberg积分公式计算加速值Tm,k m4m 1

按Romberg积分公式计算加速值

Tm,k m

4m 1T一 m 1 ,k 1 m

Tm 1,k

m=2,…，k

精度控制。对给定的精度R,若

Tm,1

Tm1,1

则终止计算，并取Tm,1为所求结果；否则返回2重复计算，直至满足要求的精度为止

程序设计：

本实验采用Matlab的M文件编写。其中待积分的函数写成 function的式，例如如下

function yy=f(x,y);

yy=x.A3;

写成如上形式即可，下面给出主程序

Romberg积分法源程序

%%cRomberg 积分法

clear

%%%分区问

b=3;

a=1;

%%%度要求

R=1e-5;

%%%用梯形公式准备初值

T(1,1)=(b-a)*(f(b)+f(a))/2;

T(1,2)=T(1,1)/2+(b-a)/2*f((b+a)/2);

T(2,1)=(4*T(1,2)-T(1,1))/(4-1);

j=2;

m=2;

%%%程序体％%%

while(abs(T(m,1)-T(m-1,1))>R); %%%度控制

j=j+1;

s=0;

for p=1:2A(j-2);

s=s+f(a+(2*p-1)*h/(2A(j-1)));

end

T(1,j)=T(1,j-1)/2+h*s/(2A(j-1)); %%% 形公式应用

for m=2:j;

k=(j-m+1);

T(m,k)=((4A(m-1))*T(m-1,k+1)-T(m-1,k))/(4A(m-1)-1);

end

%%%出Romberg积分法的函数表

I=T(m,1)

结果分析和讨论:

专业资料

进行具体的积分时，精度取 R=1e-8

100

求积分 x3dx。精确解1=

运行程序得Romberg积分法的函数表为

1.0e+007 *

4.000

3.000

2.000

由函数表知Romberg积分给出的结果为2.4999676*10^7 ,与精确没有误差, 精度很高。

求积分 131dx。精确解 I=ln3= 1.811。

运行程序得Romberg积分法的函数表为

1.333

1.667

1.068

1.303

1.846

1.559

1.111 1.000 1.535 1.048 1.027 1.130 0

1.926 1.371 1.749 1.625 1.937 0 0

1.600 1.533 1.306 1.164 0 0 0

1.313 1.850 1.179 0 0 0 0

1.593 1.046 0 0 0 0 0

1.019 0 0 0 0 0 0

从积分表中可以看出程序运行的结果为 1.019，取8位有效数字，满足要求'

求积分1四dx。

　0 x

直接按前面法进行积分，会发现系统报错，出现了。为除数的现象。出现这种情况的

原因就是当x=0时，被积函数分母出现了 0,如果用一个适当的小数（最好不要小丁程

序给定的最小误差值，但是不能小丁机器的最大精度）来代替，可以避免这个问题。本实

验取

R,可得函数表为：

0.659

0.190 0.417 0.489 0.743

0.034

0.160 0.846 0.495 0

0.368

0.092 0.138 0 0

0.722

0.726 0 0 0

0.718

0 0 0 0

故该函数的积分为0.718,取8位有效数字

4.1 2

求积分 sin x dx

本题的解析解很难给出，但运用Romberg积分可以很容易给出近似解，函数表为:

本题的解析解很难给出，但运用

0.395 0.9240.922 0.094 0.949 0.4560.1000.588 0.818 0.900 0.9590.3540.167 0.439 0.296

0.395 0.924

0.922 0.094 0.949 0.456

0.100

0.588 0.818 0.900 0.959

0.354

0.167 0.439 0.296

0.465 0.269 0.008

0.484 0.211

0 0

0.262 0

故该函数的积分为

0.262,

0 0

取8位有效数子。

0 0

结论：

Romberg积分通常要求被积函数在积分区间上没有奇点。如有奇点，且奇点为第一间断点，那么采用例3的法，还是能够求出来的，否则，必须采用其它的积分法。当然, Romberg积分的收敛速度还是比较快的。

实验报告六

题目：常微分程初值问题的数值解法

摘要：本实验主要采用经典四阶的 R-K法和四阶Adams预测-校正法来求解常微分程的

数值解。

前言：(目的和意义)

通过编写程序，进行上机计算，使得对常微分程初值问题的数值解法有更深刻的理解，掌握单步法和线性多步法是如进行实际计算的及两类法的适用围和优缺点，特别是对这两

类法中最有代表性的法：R-K法和Adams法及预测-校正法有更好的理解。通过这两种法的配合使用，掌握不同的法如配合在一起，解决实际问题。

数学原理：

对丁一阶常微分程初值问题

(Ddy/ dx f (x, y)

y(x(J y0

的数值解法是近似计算中很中的一部分。

常微分程的数值解法通常就是给出定义域上的 n个等距节点，求出所对应的函数值

yn。通常其数值解法可分为两大类：

单步法：这类法在计算yn+1的值时，只需要知道Xn+1、Xn和yn即可，就可算出。这类法典型有欧拉法和 R-K法。

多步法：这类法在计算yn+1的值时，除了需要知道xn+1、xn和扣值外，还需要知道前k步的值。典型的法如 Adams法。

f。它的具经典的R-K法是一个四阶的法。它最大的优点就是它是单步法，精度高，计算过程便丁改变步长。其缺点也很明显，计算量大，每前进一步就要计算四次函数值体的计算公式如式(

f。它的具

四阶Adams预测-校正法是一个线性多步法，它是由 Adams显式公式

yn1 yn (Ki 2K2 2K3 K4)h/6

(2)Ki f(xn,yn)

(2)

K2 f(xn h/2,yn K1h/2)

K3 f(xn h/2,yn K2h/2)

K4 f(xn h,yn K3h)

和隐式公式组成，其计算公式如式(3)所示。

预测

yn 1

yn (55 fn 59 fn 1 37 fn 2

9fn 3)h/24

(3a)

求导

fn 1

f (xn 1 , yn 1 )

(3b)

校正

yn 1

yn (9fn 1 19fn 5 fn 1 fn

2)h/24

(3c)

求导

fn 1

f g 1, yn 1)

(3d)

将局部截断误差用预测值和校正值来表示，在预测和校正的公式中分别以它们各自的阶段误差来进行弥补，可期望的到精度更高的修正的预测 -校正公式为：

预测

pn 1

(55 fn

59fn1 37 fn

2 9fn 3)h/24

修正

mn 1

p n 1

(cn

pn )251/270

求导

fn 1

f (xn

1,mn 1)

校正

cn 1

(9fn 1

19fn 5fn 1

fn 2)h/24

修正

Yn 1

cn 1

(cn 1

pn〔)19/270

求导

fn 1

f (xn

1, yn 1 )

由丁开始时无预测值和校正值可以利用，故令 po=C0=0,以后按上面进行计算。此法

的优点是可以减少计算量；缺点是它不是自开始的，需要先知道前面的四个点的值 y°,yi,y2,y3,因此不能单独使用。另外，它也不便丁改变步长。

程序设计：

本实验采用Matlab的M文件编写。其中待求的微分程写成 function的式，如下

function yy=g(x,y);

yy=-x*x-y*y;

写成如上形式即可，下面给出主程序。

经典四阶的R-K法源程序

clear

%%%长选取

h=0.1;

%%%始条件，即x=0时，y=1。

y(1)=1;

%%%解区问

a=0;

b=2;

%%%代公式

for x=a:h:b-h;

k1=g(x,y((x-a)/h+1)); %%(x) y'(x),下同。

k2=g(x+h/2,y((x-a)/h+1)+h/2*k1);

k3=g(x+h/2,y((x-a)/h+1)+h/2*k2);

k4=g(x+h,y((x-a)/h+1)+h*k3);

y((x-a)/h+2)=y((x-a)/h+1)+h*(k1+2*k2+2*k3+k4)/6;

end

四阶Adams预测-校正法源程序

%%%长选取

h=0.1;

%%%始条件

y(i)=i;

%%%解区问

a=0;

b=2;

%%%用RK迭代公式计算初始值y0,y1,y2,y3

for x=a:h:a+2*h;

k1=g(x,y((x-a)/h+1));

k2=g(x+h/2,y((x-a)/h+1)+h/2*k1);

k3=g(x+h/2,y((x-a)/h+1)+h/2*k2);

k4=g(x+h,y((x-a)/h+1)+h*k3);

y((x-a)/h+2)=y((x-a)/h+1)+h*(k1+2*k2+2*k3+k4)/6;

end

%%%用预测校正法求解

c(4)=0; %%% 正初值

p(4)=0; %%%测初值

f(1)=g(a+0*h,y(1)); g(x) y'(x),且将该值存在数组 f 中

f(2)=g(a+1*h,y (2));

f(3)=g(a+2*h,y (3));

f(4)=g(a+3*h,y (4));

for n=4:(b-a)/h;

%%%P 测

p(n+1)=y(n)+h/24*(55*f(n)-59*f(n-1)+37*f(n-2)-9*f(n-3));

%%%K 正

m(n+1)=p(n+1)+251/270*(c(n)-p(n));

%%%E 导

f(n+1)=g(a+(n+1-1)*h,m(n+1));

%%%C 正

c(n+1)=y(n)+h/24*(9*f(n+1)+19*f(n)-5*f(n-1)+f(n-2));

%%%K 正

y(n+1)=c(n+1)-19/270*(c(n+1)-p(n+1));

%%%E 导

f(n+1)=g(a+(n+1-1)*h,y(n+1));

end

结果分析和讨论：

' 2 2

计算实例一：对初值问题 y X y取步长h=0.1;计算在［-1 , 0］上的数值解。

y( 1) 0

在计算机上，运用所编的程序进行了运算，结果如下，其中第一列为在区间上的等分点，第二列为运用 R-K法的计算结果，第三列为 Adams预测-校正法计算结果。由丁本题的解析解很难求出，无法看出精度如，为此进行第二实例计算。

第一实例计算结果

-1.000 0 0

-0.000

0.240

-0.000

0.128

-0.000

0.268

-0.000

0.459

0.937

-0.000

0.451

0.732

-0.000

0.814

0.956

-0.000

0.377

0.822

-0.000

0.689

0.237

-0.000

0.114

0.911

0 0.462 0.737

计算实例二：对初值问题 y y 2X/y取步长h=0.1;计算在［0, 1.5］上的数值解。本

y(0) 1

题的解析解为y J1 2x。

同样在计算机上，运用所编的程序进行了运算，结果如下，其中第一列为在区间上的等分点，第二列为运用 R-K法的计算结果，第三列为 Adams预测-校正法计算结果，第四列为精确解。

0 1.000 1.000

1.000

0.000

1.309

1.033

0.000

1.599

1.992

0.000

1.039

1.735

0.000

1.037

1.867

1.987

0.000

1.009

1.048

1.310

0.000

1.199

1.166

1.913

0.000

1.214

1.467

1.297

0.000 1.899 1.642 1.971

0.000

1.626

1.170

1.815

1.000

1.557

1.786

1.888

1.000

1.579

1.849

1.983

1.000

1.791

1.183

1.858

1.000

1.080

1.772

1.103

1.000

1.812

1.208

1.179

1.000

2.027

1.724

2.000

根据计算结果，发现两种法的结果与精确解很接近，精度均达到 5位有效数字，但

是R-K法运算是占用的存要比 Adams预测-校正法多，即其对计算机的要求要高一些，这与理论分析吻合。当然，四阶Adams预测-校正法的启动要由R-K法给出，即景、*、％* 的值需预先知道。

相关热词搜索： 实验报告 哈工大 上机数值

哈工大-数值分析上机实验报告

最新文章

热门文章