计算方法上机实验报告.docx

计算方法上机实验报告.docx

《计算方法上机实验报告.docx》由会员分享，可在线阅读，更多相关《计算方法上机实验报告.docx（12页珍藏版）》请在冰豆网上搜索。

计算方法上机实验报告

《计算方法》上机实验报告

班级：

XXXXXX

小组成员：

XXXXXXX

XXXXXXX

XXXXXXX

XXXXXXX

任课教师：

XXX

二0—八年五月二十五日

刖言

通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton迭代法、Jacobi迭代法、

Gauss-Seidel迭代法、Newton插值法、Lagrange插值法和Gauss

求积公式等六种算法的原理和使用方法，并参考课本例题进行了

MATLAB程序的编写。

以下为本次上机实验报告，按照实验内容共分为六部分。

实验一：

一、实验名称及题目:

Newton迭代法

X3IXI1

例2.7（P38）应用Newton迭代法求

在-=1附近的数值

解:

并使其满足二、解题思路：

设X'是f（x）0的根选取X。

作为X'初始近似值，过点Xo,f（Xo）做曲

线yf（x）的切线L丄的方程为yf（Xo）f'（x°）（xx。

）,求出L与x轴交

点的横坐标X1xo鴿称X1为X'的一次近似值过点风心））做曲线

yf（x）的切线，求该切线与X轴的横坐标X2Xi卫^称X2为X'的二次f'（Xi）

近似值重复以上过程，得X'的近似值序列Xn,把XniXn竺称为x'

f'（Xn）

的n1次近似值，这种求解方法就是牛顿迭代法。

三、Matlab程序代码:

functionnewton_iteration（xO,tol）symsz%定义自变量formatlong%定义精度f=z*z*z-z-1;

f1=diff（f）;%求导y=subs（f,z,x0）;

y1=subs（f1,z,x0）;%向函数中代值x1=x0-y/y1;k=1;

whileabs（x1-x0）>=tol

x0=x1;

y=subs（f,z,x0）;y1=subs（f1,z,x0）;x1=x0-y/y1;k=k+1;

end

x=double（x1）

K

四、运行结果：

>>tlevfn-n_itezatinnll.O.00000001）

x=

1.»47L7»B7244746

实验二：

一、实验名称及题目:

Jacobi迭代法

例3.7（P74）:

试利用Jacobi迭代公式求解方程组

要求数值解

为方程组的精确解.

二、解题思路：

首先将方程组中的系数矩阵A分解成三部分，即：

ALDU，D为对角阵，L为下三角矩阵，U为上三角矩阵。

之后确定迭代格式,X（k1）b*x（k）f，（k0,1,2,k即迭代次数），B称为迭代矩阵。

最后选取初始迭代向量X（0），开始逐次迭代。

最后验证精度。

（迭代阵：

X（k1）DUX（k）D1b。

）

雅克比迭代法的优点明显，计算公式简单，每迭代一次只需计算一次矩阵和向量的乘法，且计算过程中原始矩阵A始终不变，比较容易并行计算。

然而这种迭代方式收敛速度较慢，而且占据的存储空间较大。

三、Matlab程序代码:

functionjacobi（A,b,xO,eps,x1）

D=diag（diag（A））;%求A的对角矩阵

L=-tril（A,-1）;%求A的下三角矩阵

U=-triu（A,1）;%求A的上三角矩阵

B=D\（L+U）;

f=D\b;

x=B*xO+f;

n=1;%迭代次数

whilenorm（x-x1）>=eps

x=B*x+f;

n=n+1;

end

formatlong

n

x

jingdu=norm（x-x1）

四、运行结果：

»』-1-I-L-l14-I-l.-l-1-1WWIS-4,[I23

*

实验三:

一、实验名称及题目:

Gauss-Seidel迭代法

例3.8（P75）:

试利用Gauss-Seidel迭代公式求解方程组

并使其

数值解

为方程组的精确解

解题思路:

Gauss-Seidel迭代法与Jacobi迭代法思路相近，首先将方程组中的系数矩阵A分解成三部分，即：

ALDU，D为对角阵，L为下三角矩阵，U为上三角矩阵。

之后确定迭代格式，X（k1）B*X（k）f，（k0,1,2,k即迭代次数），B称为迭代矩阵。

最后选取初始迭代

向量X0，开始逐次迭代。

最后验证精度。

（迭代阵：

X（k1）（DL）1UX（k）（DL）1b。

）

Gauss-Seidel迭代法与Jacobi迭代法相比速度更快，但不全如此。

有例子表明：

Gauss-Seidel迭代法收敛时，Jacobi迭代法可能不收敛；而Jacobi迭代法收敛时，Gauss-Seidel迭代法也可能不收敛。

三、Matlab程序代码:

functiongauss_seidel（A,b,xO,eps,x1）D=diag（diag（A））;%求A的对角矩阵

L=-tril（A,-1）;%求A的下三角矩阵

U=-triu（A,1）;%求A的上三角矩阵

B=（D-L）\U;

f=（D-L）\b;

x=B*xO+f;

n=1;%迭代次数

whilenorm（x1-x）>=eps

x=B*x+f;n=n+1;

end

formatlong

n

x

jingdu=norm（x1-x）

四、运行结果:

»■ME.sndflKCB-1-3-l.-l10-1-1;-1-1B-1-1-1-110]R[-4.12.8.34],[O.OjO.OL10^-€[1.3j;4]）n-

8

x■

0.躬曲討；E刃灯2開

』・M9^BlB1194SL0]i

2.H9«7>2eiW3U

jincdu=

实验四:

一、实验名称及题目:

Lagrange插值法

例4.1（P88）:

给定函数I--丄沁及插值节点

_71_7L_37T_71

'°-1J''一斗"三一Jr_丄试构造Lagrange插值多项式，

给出其误差估计，并由此计算「匸）及其误差.

二、解题思路：

一般来说，如果我们有n个点x1,y1,...,xn,yn，各人互不相同。

那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：

n

L（x）yjlj（x），其中每个lj（x）为拉格朗日基本多项式（或称插值基函

j0

数），其表达式为:

lj（x）

Xi

。

i0,ijXj

Xi

三、Matlab程序代码:

functiony=lagrange（xO,x）n=length（xO）;%向量长度s=0;

fork=1:

n%k从1至Un的循环p=1.0;forj=1:

n

ifj〜=k%~=”不等于的意思

p=p*（x-x0（j））/（x0（k）-x0（j））;

end

endy0=x0（k）*（1+cos（x0（k）））;s=p*y0+s;

end

formatlong

s

wucha=abs（x*（1+cos（x））-s）

四、运行结果:

»lagr«ngt（[0pi/8pi/43*pi/8pi/2]»3*pi/16）

£=

.0791********

wucha-

3.33S58S31BL7o972«-04

五、Lagrange插值图像绘制

%Lagrange插值图像算法

x=linspace（0,1002,200）;

s=linspace（0,1000,200）;

x0=[0;pi/8;pi/4;3*pi/8;pi/2];

n=length（x0）;

s=0;

fork=1:

n

p=1.0;

forj=1:

n

ifj~=k

p=p.*（x-x0（j））/（x0（k）-x0（j））;

end

end

y0=x0（k）*（1+cos（x0（k）））;

s=p*y0+s;

end

plot（x,s,'r'）;

gridon;

title（'Lagrange2?

?

ui?

?

?

）

xlabel（'X'）,ylabel（'Y'）;

axisnormal;

实验五:

一、实验名称及题目:

Newton插值法

例4.3（P96）:

已知，试取插值节点

Xq=0.35j=O.SOj%2=0=65』丸寻崔=0<95构造4次

Newton插值多项式，由此计算W亦的逼近值，并指出其绝对误差.

二、解题思路：

nnn

将拉格朗日插值公式Ln（x）三4yjlj（x）yj中的Ln（x）改

写成：

Nn（x）a。

ai（xx°）

j0i0XjXij0

xxf（x）f（xXo）

x,x0?

xXo

a2（Xx0）（Xx1）...an（xxo）...（xxnl）,其中，

a。

®,an为待定定系数。

又a。

f（x°）,ai

anfX,Xo,...,XnfXN’...,xn1fx°,心…,x.。

彳将a°,ai,an带入Nn（x）可

xXn

得：

f（x）f（Xo）f[Xo,Xi]（XXo）...f[Xo,Xi,...,Xn]（XX°）（X为）...（XXn1）。

三、Matlab程序代码：

functionnewton_interpolation（xO,x）formatlongn=length（xO）;

symszf=sqrt（1+cosh（z）A2）;a

（1）=subs（f,z,x0

（1））;fork=1:

n-1

yO=subs（f,z,xO（k））;

y仁subs（f,z,x0（k+1））;

d（k,1）=（y1-yO）/（xO（k+1）-xO（k））;k阶差商

end

forj=2:

n-1

%二阶

fork=1:

n-jd（k,j）=（d（k+1,j-1）-d（k,j-1））/（x0（k+j）-x0（k））;差商及以上

end

end

double（d）

forj=2:

n

a（j）=d（1,j-1）;

end

b

（1）=1;c

（1）=a

（1）;

forj=2:

n

b（j）=（x-x0（j-1））.*b（j-1）;c（j）=a（j）.*b（j）;

end

np=double（sum（c））wucha=double（abs（np-subs（f,z,x）））

四、运行结果：

>Jini4rpAlBllifrnt39.0.SQLAL49.0.HOL0.0・T）

4亠Ou44£B2aTaU03tli

40.34LZftII322«0»

D-i5Zl!

!

7r€iT4SEE2®^5O.G21

Qu“

thia«TM«i3LlS!

l443aHTH003>3«2«S3

0Q

00

■

1.M4I319T13173747

■rucHa■

IBI«S42S»-07

五、Newton插值图像绘制

实验六：

一、实验名称及题目:

Gauss求积公式

例5.7（P140）:

试构造Gauss型求积公式

At

（I+O3

fMdxIylo/Uo）+血f伽）+如fg）

并由此计算积分

解题思路:

AJ0

A2f

1

设高斯-勒让德求积公式是：

if（x）dxAgf

-1

依次代入f（x）1,x,x2,解得A0A29，A8。

利用换元公式

宁x专dx变换原式的积分上下限’在套用高斯

-勒让德求积公式求得积分

三、Matlab程序代码:

functiongauss（a,b）

symstf=sqrt（t）/（1+t）A2;

P=[-sqrt（3/5）0sqrt（3/5）];

A=[5/98/95/9];

s=0;

fori=1:

3;

x=P（i）;

y=subs（f,t,（b-a）*x/2+（a+b）/2）;s=s+A（i）*y;

end

formatlong

S=double（（a-b）/2*s）

四、运行结果：

1-

、、iauss<1,0）

结束语

在本学期的《计算方法》课程学习中，我们感受到了巧妙的数学计算方法在解决实际问题中带来的便利与高效，借助计算机解决科学计算问题也是我们当代大学生应当掌握的必要技能。

在本课程的上机实验过程中，我们亲自体验了课本中所学到的算法在计算机上的使用，使用计算机语言，如MATLAB和C语言等进

行编程，加深了我们对各类算法及其理论的理解，并进一步激发了我们对《计算方法》这门课程进行持续学习的学习兴趣