平行平面腔自再现模的matlab数值计算

平行平面腔自再现模的数值解法

石鹏（2010302749）

目前，平行平面腔在中等以上激光器中仍普遍应用，但是平行平面腔的自再现模积分方程至今尚未得到精确的解析解。本文就条形腔和矩形对积分方程做数值计算。计算的思想基于Fox-Li 数值迭代法，计算用matlab 软件实现，结果可以很直观地用图形表示出来。

1.矩形腔的分离变量法

先来看谐振腔的自再现模。所谓自再现，就是指不管初始分布1u 的具体特征如何，只要经过足够多次的往返渡越后，所生成的场分布都将明显带有衍射的痕迹，而且，具有这种特点的场分布将不再受到衍射作用的影响，形成一种稳定的场分布。

设(),q u x y ''为经过q 次渡越后在某一镜面上的光场分布，根据菲涅耳—基尔霍夫衍射积分公式，可以得到光波再渡越一次腔长距离后的光场分布1q u +与q u 之间的关系为：

()()

()1,,1cos 4ik q q

ik

e u x y u x y ds ρ

θπ

ρ

-+'''=

+?? （1）

按照自再现模的概念，1q u +和q u 应满足关系：1q q u u σ+=。光场分布q u 的模长为振幅分布，相角为相位分布。

对于矩形腔，ρ的计算公式为：

ρ=

按幂级数展开近似为：()()2

2

22x x y y L L

L

ρ''--=+

+

代入（1），可以得到积分方程的核为：()()()2222,,,x x y y ik L L ikL i

K x y x y e e

L

λ??

''--??+

??-??

''=

可以看出，方形镜平行平面腔的自再现模的积分方程对x 、y 两个坐标是对称的。令

()()(),mn m n u x y u x u y =，mn m n σσσ= （2）

这样就可以对积分方程进行分离变量，x 和y

的方程具有相同的形式：

()()()2

2x x a ik

L

m m m a u x e

u x dx σ'---''= （3）

这就是条形腔的积分方程，求解矩形腔的自再现模的关键就是解这个方程。

2.积分方程的数值解法

积分方程（1），（3）在数学上称为第二类弗里德霍姆方程。这儿使用Fox-Li 数值迭

代的方法来求解。方程（3）可以表示为以下形式：

()()()1a

q q a

u x K x x u x dx +-'''=-? （4）

(

)()2

2x x ik

L

K x x '--'-=

（5）

其中()K x x '-为积分的核函数。方程（4）的含义是场分布()q u x 经过一次传输以后得到新的场分布()1q u x +。根据自再现模的性质，不管初始分布()1u x 具有什么样的形式，用方程（4）经过足够多次迭代得到的光场分布都具有相同的形式。

现在任意给一个初始分布，来解方程（4）。先将方程（4）离散化：

()()()

11

2n

q i q i i a u x K x x u x n +=''=-∑ （6） 再次迭代就要用到()1q i u x +'，把i x '代入上式，令i 分别等于1,2,3，计算()

1q i u x +'的前

几项，()()()11112n q i q i i a u x K x x u x n +=''=-∑；()()()

1221

2n

q i q i i a

u x K x x u x n +=''=-∑；()()()

1331

2n q i q i i a u x K x x u x n +=''=-∑，可以发现（6）式具有矩阵的形式，把它写出来： ()()()()()

()()()()()()()()()()11111121122212221n 122q q n q q n q n q n n n n u x u x K x x K x x K x x u x u x K x x K x x K x x a n u x u x K x x K x x K x x +++????

---?? ?

? ?--- ? ?

?= ? ? ? ? ? ? ? ?

?---??????

……………………

(7)

把上式简化表示为1q q U KU +=，可以发现矩阵表示有特殊的含义：在（4）中，()K x x '-是一个关于x 和x '的函数，而在（7）中K 是一个数值矩阵。传输一次就是给q U 左乘一个K ，于是，初始分布经过q 次传输以后就是给初始分布左乘一个q

K ：

1q q U K U = （8）

这就是积分方程（4）的数值方法的迭代公式。可以发现，用积分方程迭代的重点就是

计算核矩阵K 。我把它做成了一个函数文件kernel ：

function K = kernel( Lambda , L , x )

%KERNEL 计算积分方程的核K

k = 2*pi/Lambda ;

[ X , Y ] = meshgrid(x) ;

K = sqrt( 1i/(Lambda*L) * exp(-1i*k*L) ) * exp(-1i*k * (Y-X).^2 / (2*L) ) ;

end

Kernel函数的输入参数为波长、腔长和一个把坐标轴离散化后的向量x 。其中用到了一个生成网格数据的函数meshgrid，它可以得到一个横向向右递增的

矩阵X和一个纵向向下递增的矩阵Y，用来计算矩阵K。这种算法的结果得到了

一个数值矩阵，并且把函数积分的N次迭代转化为对数值矩阵求N次方。matlab 计算矩阵的速度比计算函数的速度要快，这就增加了程序的效率，优化了算法。

这儿得到的K中并没有步长2a/n，因为n是在函数外部定义的，所以在函数外部调用比较方便。

3.条形腔的自再现模

方程（3）就是条形腔的自再现模的积分方程，经过前面的讨论，它的解法就是利用公式（8）来进行迭代。计算条形腔的振幅和相位分布的函数为bar1，下面给出bar1函数的一段代码：

M = 500 ; T = 300 ;

x = linspace(-a,a,M) ; x1 = ones( M , 1 ) ;

K = kernel( Lambda , L , x ) ;

y = (2*a/(M-1) * K)^T *x1 ;

y0 = max(abs(y)) ;

y1 = abs(y)/y0 ; y2a = angle(y) ;

y2 = xiugaiAngle( y2a , M ) ;

其中，M为离散程度，T为传播次数，即迭代次数。代码中有一个全一列向量x1，这就是我给的初始分布，全一表示它是一个平面波。第三行生成了一个需要的K矩阵，真正进行迭代计算的是第四行。步长2*a/(M-1)在这儿调用，迭代了T次。这儿输出的y就是最终的光场分布。最后得到的y1为振幅分布，y2为相位分布。

在计算相位分布时没有直接输出y2a ，而是用到了一个函数xiugaiAngle ，可以称它为修改相位函数。因为我们发现在计算相位时，angle 函数的值域为()-ππ，，这就意味着当相位超过π

时angle 函数会自动跳到-π，这样的话画出来的相位分布就会有一个个“坑”，如右图：

要做的修改就是找出这些“坑”并且给它加上2π 。xiugaiAngle 函数如下：

function y = xiugaiAngle( x , M ) y = x ;

x1 = [ x(1) ; x ] ; x2 = [ x ; x(M) ] ; xa = x2 - x1 ;

y1 = find(xa < -2*pi + 0.2 ) ; y2 = find( xa > 2*pi - 0.2 ) ; n = length( y1 ) ; for ii = 1 : n

y = [ y(1:y1(ii)-1) ; y(y1(ii):y2(ii)-1)+2*pi ; y(y2(ii):M) ] ; end end

用find 函数找到“坑”的位置，用for 循环给它加上2π 。

4.矩形腔的自再现模

从方程（1）由（2）到方程（3），说明了矩形腔的分离变量法。现在解出了方程（3），那么就可以由方程（2）到方程（1），计算矩形腔的自再现模。计算用的函数是rectangleq2 。

下面一段代码是计算矩形腔的振幅分布的：

Z = xa*conj(ya)' ;

Z1 = abs(Z)/max( max(abs(Z)) ) ;

其中xa和ya分别代表了x方向的光场分布和y方向上的光场分布。它们都是列向量，对ya转置（matlab中的算符‘ 表示取一个矩阵的共轭转置，还要再求一次共轭才能得到转置），代码中的乘法是矩阵乘法，也就是说乘得的结果是一个矩阵，也就是由公式（2）得到的光场分布。

由于xiugaiAngle函数是一维的，所以不能用上面得到的Z来计算相角分布。相角分布的计算由下面一段代码实现：

xa2 = angle(xa) ; ya2 = angle(ya) ;

xb2 = xiugaiAngle( xa2 , M ) ; yb2 = xiugaiAngle( ya2 , M ) ;

[ X2 , Y2 ] = meshgrid( xb2 , yb2 ) ;

Z2 = X2 + Y2 ;

复数求积后相角是相加的，上面第四行代码就是实现相角相加。

参考文献：

[1]李东激光原理与技术课件

[2]俞宽新激光原理与机关技术北京工业大学出版社149-164

[3]孙蓬matlab基础教程

[4]吴鹏MATLAB高效编程技巧与应用：25个案例分析北京航空航天大学出版社98-99

说明：一共有7个m文件：

bar000.m 脚本文件是这道题的第一个成功的m文件，画的是条形腔的振幅分布

d = 1 ; k = 2*pi/d ; a = 25*d ; L = 100*d ;

M = 300 ; T = 300 ;

x = linspace(-a,a,M) ; x1 = ones( M , 1 ) ;

[ X , Y ] = meshgrid(x) ;

K1 = sqrt( 1i/(d*L) * exp(-1i*k*L) ) ;

K = K1 * exp(-1i*k * (Y-X).^2 / (2*L) ) ;

y1 = (2*a/(M-1) * K)^T *x1 ;

y = abs(y1)/max(abs(y1)) ;

plot(x,y)

bar1.m 函数文件用来计算并绘制条形腔的振幅和相位分布function bar1( Lambda , a , L )

%Bar1 绘制条形腔的振幅分布和相位分布

% 输入参数为波长，镜面宽度，腔长

M = 500 ; T = 300 ; %M为离散程度，T为传播次数x = linspace(-a,a,M) ; x1 = ones( M , 1 ) ;

K = kernel( Lambda , L , x ) ;

y = (2*a/(M-1) * K)^T *x1 ;

y0 = max(abs(y)) ;

y1 = abs(y)/y0 ; y2a = angle(y) ;

y2 = xiugaiAngle( y2a , M ) ;

na_L = int2str(L) ;

na_a = int2str(a) ;

na_M = int2str(M) ; na_T = int2str(T) ;

namez = [ '条形腔振幅分布' ] ;

namex = [ '条形腔相位分布' ] ;

name_xq = [ 'a：' na_a '倍波长' 'L：' na_L '倍波长' ] ;

name_xh = [ '离散程度：' na_M '迭代次数：' na_T ] ;

figure(1) ; plot(x,y1)

title( [ namez char(10) name_xq char(10) name_xh] ) ;

saveas(gcf,[ namez na_a '_' na_L '.jpg' ])

figure(2) ; plot(x,y2)

title( [ namex char(10) name_xq char(10) name_xh] ) ;

saveas(gcf,[ namex na_a '_' na_L '.jpg' ])

end

barq.m 脚本文件调用bar1.m来实现条形腔光场的计算

disp('绘制条形腔的振幅分布和相位分布') ;

Lambda = input('请输入波长λ= ') ;

a = input('请输入镜面宽度2a = ') /2 ;

L = input('请输入腔长L = ') ;

bar1( Lambda , a , L ) ;

kernel.m 函数文件计算积分核K

function K = kernel( Lambda , L , x )

%KERNEL 计算积分方程的核K

% 可以输入参数为波长和腔长

k = 2*pi/Lambda ;

[ X , Y ] = meshgrid(x) ;

K = sqrt( 1i/(Lambda*L) * exp(-1i*k*L) ) * exp(-1i*k * (Y-X).^2 / (2*L) ) ; End

rectangleq.m 脚本文件调用rectangle2.m计算矩形腔的光场分布disp('绘制矩形腔的振幅分布和相位分布') ;

Lambda = input('请输入波长λ= ') ;

a = input('请输入镜面长度2a = ') /2 ;

b = input('请输入镜面宽度2b = ') /2 ;

L = input('请输入腔长L = ') ;

rectangleq2( Lambda , a , b , L )

rectangleq2.m 函数文件计算矩形腔的光场分布

function rectangleq2( Lambda , a1 , b1 , L1 )

%RECTANGLE2 绘制矩形腔的振幅分布和相位分布

% 修补相位分布的缺陷

a = a1/Lambda ;

b = b1/Lambda ; L = L1/Lambda ;

M = 500 ; T = 300 ;

x = linspace(-a,a,M) ; y = linspace(-b,b,M) ;

x0 = ones( M , 1 ) ;

[ X , Y ] = meshgrid( x , y ) ;

Kx = kernel( Lambda , L , x ) ;

xa = (2*a/(M-1) * Kx)^T *x0 ;

Ky = kernel( Lambda , L , y ) ;

ya = (2*b/(M-1) * Ky)^T *x0 ;

Z = xa*conj(ya)' ;

Z1 = abs(Z)/max( max(abs(Z)) ) ;

%计算相位分布

xa2 = angle(xa) ; ya2 = angle(ya) ;

xb2 = xiugaiAngle( xa2 , M ) ; yb2 = xiugaiAngle( ya2 , M ) ;

[ X2 , Y2 ] = meshgrid( xb2 , yb2 ) ;

Z2 = X2 + Y2 ;

na_L = int2str(L) ;

na_a = int2str(a) ; na_b = int2str(b) ;

na_M = int2str(M) ; na_T = int2str(T) ;

namez = '矩形腔振幅分布' ;

namex = '矩形腔相位分布' ;

name_xq = [ 'a：' na_a '倍波长' ' b：' na_b '倍波长' 'L：' na_L '倍波长' ] ; name_xh = [ '离散程度：' na_M '迭代次数：' na_T ] ;

figure(1) ; mesh( X , Y , Z1 )

title( [ namez char(10) name_xq char(10) name_xh] ) ;

saveas(gcf,[ namez na_a '×' na_b '.jpg' ])

figure(2) ; mesh( X , Y , Z2 )

title( [ namex char(10) name_xq char(10) name_xh] ) ;

saveas(gcf,[ namex na_a '×' na_b '.jpg' ])

end

xiugaiAngle.m 函数文件修改相位函数

function y = xiugaiAngle( x , M )

%UNTITLED2 Summary of this function goes here

% Detailed explanation goes here

y = x ;

x1 = [ x(1) ; x ] ; x2 = [ x ; x(M) ] ;

xa = x2 - x1 ;

y1 = find(xa < -2*pi + 0.2 ) ;

y2 = find( xa > 2*pi - 0.2 ) ;

n = length( y1 ) ;

for ii = 1 : n

y = [ y(1:y1(ii)-1) ; y(y1(ii):y2(ii)-1)+2*pi ; y(y2(ii):M) ] ;

end

end

有11张图片，图片名就是图片的说明

数值计算方法实验指导(Matlab版)

《数值计算方法》实验指导 （Matlab 版） 肇庆学院数学与统计学学院 计算方法课程组

1. 实验名称 实验1 算法设计原则验证(之相近数相减、大数吃小数和简化计算步骤) 2. 实验题目 有效数字的损失. 123 )与1000个较小的数(3 10 15)的和，验证 大数吃小数的现象. (3)分别用直接法和秦九韶算法计算多项式 P(x) a 0x n a 1x n 1 在x =1.00037 处的值?验证简化计算步骤能减少运算时间. n 1 对于第(3)题中的多项式P (x )，直接逐项计算需要n (n 1) 2 1 次乘法 和n 次加法，使用秦九韶算法 P(x) (((a °x ajx a 2)x a . 则只需要n 次乘法和n 次加法. 3. 实验目的 验证数值算法需遵循的若干规则. 4. 基础理论 设计数值算法时，应避免两个相近的数相减、防止大数吃小数、简化计算步骤减少运算 次数以减少运算时间并降低舍入误差的积累. 两相近的数相减会损失有效数字的个数， 用一 《数值计算方法》实验 1报告 班级： 20xx 级 XXXXx 班 学号： 20xx2409xxxx 姓名： XXX 成绩: ⑴取 z 1016，计算z 1 Z 和 1/（、z 1 Z)，验证两个相近的数相减会造成 (2)按不同顺序求一个较大的数( a n 1 X a n

个大数依次加小数，小数会被大数吃掉，乘法运算次数太多会增加运算时间. 5.实验环境 操作系统：Win dows xp ；程序设计语言：Matlab 6.实验过程 (1)直接计算并比较； (2)法1 :大数逐个加1000个小数，法2 :先把1000个小数相加再与大数加； (3)将由高次项到低次项的系数保存到数组A[n]中，其中n为多项式次数. 7.结果与分析 (1)计算的~1V Z = _______________________________ ,1/( ~1 < z) ____________________ . 分析： (2)123逐次加1000个3 10 6的和是_________________________ ，先将1000个3 10 6相 加，再用这个和与123相加得_______________________ . 分析： (3)计算__________ 次的多项式： 直接计算的结果是___________________ ，用时___________________ ; 用秦九韶算法计算的结果是____________________ ，用时 ___________________ 分析：

数值分析Matlab作业

数值分析编程作业

2012年12月 第二章 14.考虑梯形电阻电路的设计，电路如下： 电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组： 12 123 234 345 456 567 678 78 22/ 2520 2520 2520 2520 2520 2520 250 i i V R i i i i i i i i i i i i i i i i i i i i -= -+-= -+-= -+-= -+-= -+-= -+-= -+= 这是一个三对角方程组。设V=220V，R=27Ω，运用追赶法，求各段电路的电流量。Matlab程序如下： function chase () %追赶法求梯形电路中各段的电流量 a=input('请输入下主对角线向量a='); b=input('请输入主对角线向量b='); c=input('请输入上主对角线向量c='); d=input('请输入右端向量d='); n=input('请输入系数矩阵维数n='); u(1)=b(1); for i=2:n l(i)=a(i)/u(i-1); u(i)=b(i)-c(i-1)*l(i); end y(1)=d(1); for i=2:n y(i)=d(i)-l(i)*y(i-1); end x(n)=y(n)/u(n); i=n-1; while i>0 x(i)=(y(i)-c(i)*x(i+1))/u(i); i=i-1; end x 输入如下:

请输入下主对角线向量a=[0,-2,-2,-2,-2,-2,-2,-2]; 请输入主对角线向量b=[2,5,5,5,5,5,5,5]; 请输入上主对角线向量c=[-2,-2,-2,-2,-2,-2,-2,0]; 请输入方程组右端向量d=[220/27,0,0,0,0,0,0,0]; 请输入系数矩阵阶数n=8 运行结果如下： x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477 第三章 14.试分别用(1)Jacobi 迭代法；(2)Gauss-Seidel 迭代法解线性方程组 1234510123412191232721735143231211743511512x x x x x ?????? ??????---????????????=--?????? --?????? ??????---?????? 迭代初始向量 (0)(0,0,0,0,0)T x =。 （1）雅可比迭代法程序如下： function jacobi() %Jacobi 迭代法 a=input('请输入系数矩阵a='); b=input('请输入右端向量b='); x0=input('请输入初始向量x0='); n=input('请输入系数矩阵阶数n='); er=input('请输入允许误差er='); N=input('请输入最大迭代次数N='); for i=1:n for j=1:n if i==j d(i,j)=a(i,j); else d(i,j)=0; end end end m=eye(5)-d\a; %迭代矩阵 g=d\b; x=m*x0+g; k=1; while k<=N %进行迭代 for i=1:5 if max(abs(x(i)-x0(i))) >er x=m*x+g; k=k+1;

数值分析的matlab实现

第2章牛顿插值法实现 参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55. 求牛顿插值多项式和差商的MA TLAB 主程序： function[A,C,L,wcgs,Cw]=newpoly(X,Y) n=length(X);A=zeros(n,n);A(:,1) =Y'; s=0.0;p=1.0;q=1.0;c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1; end C=A(n,n);b=poly(X(n));q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C);Q=poly2sym(q1); syms M wcgs=M*Q/c1;Cw=q1/c1; （1）保存名为newpoly.m 的M 文件 （2）输入MA TLAB 程序 >> X=[242,243,249,250]; >> Y=[13.681,13.526,13.098,13.095]; >> [A,C,L,wcgs,Cw]=newpoly(X,Y) 输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其 ) ()()1(ξ+n n f x R 的系数向量Cw 。 A = 13.6810 0 0 0 13.5260 -0.1550 0 0 13.0980 -0.0713 0.0120 0 13.0950 -0.0030 0.0098 -0.0003 C = 1.0e+003 *

matlab实现数值计算功能源程序(个人整理)

matlab数值计算功能 1,基础运算 （1）多项式的创建与表达 将多项式(x-6)(x-3)(x-8)表示为系数形式 a=[6 3 8] % 写成根矢量 pa=poly(a)% 求出系数矢量 ppa=poly2sym(pa,'x') % 表示成符号形式 ezplot(ppa,[-50,50]) 求3介方阵A的特征多项式 a=[6 2 4;7 5 6;1 3 6 ]; pa=poly(a)% 写出系数矢量 ppa=poly2sym(pa) %表示成符号形式 ezplot(ppa,[-50,50]) % 绘图 求x^3-6x^2-72x-27的根。 a=[1,-6,-72,-85]; % 写出多项式系数矢量 r=roots(a) % 求多项式的根 （2）多项式的乘除运算 c=conv(a,b) %乘法 [q,r]=deconv(c,a)% 除法 求a(s)=s^2+2s+3乘以b(s)=4s^2+5s+6的乘积 a=[1 2 3] b=[4 5 6] % 写出系数矢量 c=conv(a,b) c=poly2sym(c,'s') % 写成符号形式的多项式 展开(s^2+2s+2)(s+4)(s+1)并验证结果被(s+4),(s+3)除后的结果。c=conv([1,2,2],conv([1,4],[1,1])) cs=poly2sym(c,'s') c=[1 7 16 18 8] [q1,r1]=deconv(c,[1,4]) [q2,r2]=deconv(c,[1,3]) cc=conv(q2,[1,3]) test=((c-r2)==cc)

其他常用的多项式运算命令 pa=polyval(p,s) % 按数组规则计算给定s时多项式的值 pm=polyvalm(p,s)% 按矩阵规则计算给定s时多项式的值 [r,p,k]=residue(b,a) % 部分分式展开，b,a分别是分子，分母多项式系数矢量。r,p,k分别是留数，极点和值项矢量。 p=poly(x,y,n) % 用n介多项式拟合给定的数据 polyder(p) %多项式微分 x=[1 2 3 4 5]; y=[5.5 43.1 128 290.7 498.4]; p=polyfit(x,y,3) x2=1:0.1:5; y2=polyval(p,x2); plot(x,y,'o',x2,y2) 2，线性代数 1，求解方程的根 t=[0 0.3 0.8 1.1 1.6 2.3]'; y=[0.82 0.72 0.63 0.60 0.55 0.50]'; e=[ones(size(t)) exp(-t)] c=e\y t1=[0:0.1:2.5]'; y1=[ones(size(t1)),exp(-t1)]*c; plot(t1,y1,'b',t,y,'ro') 2，逆矩阵及行列式 inv(a) det(a) pinv(a) 3，矩阵分解（略） 4，数据分析 max(x)求x各列的最大元素 mean(x)求x各列的平均值 median(x)找出x各列的中位元素 min(x)求出x各列的最小元素 S=cumsum(x)求x各列元素累计和 sort(x)使x的各列元素按递增排序。 std(x)求x各列的标准差。 sum(x)求x各列元素之和 prod(x)求x各列元素之积

第2讲 matlab的数值分析

第二讲MATLAB的数值分析 2-1矩阵运算与数组运算 矩阵运算和数组运算是MATLAB数值运算的两大类型，矩阵运算是按矩阵的运算规则进行的，而数组运算则是按数组元素逐一进行的。因此，在进行某些运算（如乘、除）时，矩阵运算和数组运算有着较大的差别。在MATLAB中，可以对矩阵进行数组运算，这时是把矩阵视为数组，运算按数组的运算规则。也可以对数组进行矩阵运算，这时是把数组视为矩阵，运算按矩阵的运算规则进行。 1、矩阵加减与数组加减 矩阵加减与数组加减运算效果一致，运算符也相同，可分为两种情况： （1）若参与运算的两矩阵（数组）的维数相同，则加减运算的结果是将两矩阵的对应元素进行加减，如 A=[1 1 1;2 2 2;3 3 3]; B=A; A+B ans= 2 2 2 4 4 4 6 6 6 （2）若参与运算的两矩阵之一为标量（1*1的矩阵），则加减运算的结果是将矩阵（数组）的每一元素与该标量逐一相加减，如 A=[1 1 1;2 2 2;3 3 3]; A+2 ans= 3 3 3 4 4 4 5 5 5 2、矩阵乘与数组乘 （1）矩阵乘 矩阵乘与数组乘有着较大差别，运算结果也完全不同。矩阵乘的运算符为“*”，运算是按矩阵的乘法规则进行，即参与乘运算的两矩阵的内维必须相同。设A、B为参与乘运算的 =A m×k B k×n。因此，参与运两矩阵，C为A和B的矩阵乘的结果，则它们必须满足关系C m ×n 算的两矩阵的顺序不能任意调换，因为A*B和B*A计算结果很可能是完全不一样的。如：A=[1 1 1;2 2 2;3 3 3]; B=A;

A*B ans= 6 6 6 12 12 12 18 18 18 F=ones(1,3); G=ones(3,1); F*G ans 3 G*F ans= 1 1 1 1 1 1 1 1 1 （2）数组乘 数组乘的运算符为“.*”，运算符中的点号不能遗漏，也不能随意加空格符。参加数组乘运算的两数组的大小必须相等（即同维数组）。数组乘的结果是将两同维数组（矩阵）的对应元素逐一相乘，因此，A.*B和B.*A的计算结果是完全相同的，如： A=[1 1 1 1 1;2 2 2 2 2;3 3 3 3 3]; B=A; A.*B ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 B.*A ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 由于矩阵运算和数组运算的差异，能进行数组乘运算的两矩阵，不一定能进行矩阵乘运算。如 A=ones(1,3); B=A; A.*B ans= 1 1 1 A*A ???Error using= =>

数值分析的MATLAB程序

列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1);

end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x')

同济大学数值分析matlab编程题汇编

MATLAB 编程题库 1.下面的数据表近似地满足函数2 1cx b ax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像. 625 .0718.0801.0823.0802.0687.0606.0356.0995 .0628.0544.0008.0213.0362.0586.0931.0i i y x ---- 解： x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y; a=z(1); b=z(2); c=z(3); xh=-1:0.1:1; yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*')

2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数 10 的近似根，并写出调用方式： 精度为10 解： >> edit gexianfa.m function [x iter]=gexianfa(f,x0,x1,tol) iter=0; while(norm(x1-x0)>tol) iter=iter+1; x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end >> edit f.m function v=f(x) v=x.*log(x)-1; >> edit g.m function z=g(y) z=y.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8

第06章_MATLAB数值计算_例题源程序汇总

第6章 MATLAB 数值计算 例6.1 求矩阵A 的每行及每列的最大和最小元素，并求整个矩阵的最大和最小元素。 1356 78256323578255631 01-???? -? ?=???? -??A A=[13,-56,78;25,63,-235;78,25,563;1,0,-1]; max(A,[],2) %求每行最大元素 min(A,[],2) %求每行最小元素 max(A) %求每列最大元素 min(A) %求每列最小元素 max(max(A)) %求整个矩阵的最大元素。也可使用命令：max(A(:)) min(min(A)) %求整个矩阵的最小元素。也可使用命令：min(A(:)) 例6.2 求矩阵A 的每行元素的乘积和全部元素的乘积。 A=[1,2,3,4;5,6,7,8;9,10,11,12]; S=prod(A,2) prod(S) %求A 的全部元素的乘积。也可以使用命令prod(A(:)) 例6.3 求向量X =(1！，2！，3！，…，10！)。 X=cumprod(1:10) 例6.4 对二维矩阵x ，从不同维方向求出其标准方差。 x=[4,5,6;1,4,8] %产生一个二维矩阵x y1=std(x,0,1) y2=std(x,1,1) y3=std(x,0,2) y4=std(x,1,2) 例6.5 生成满足正态分布的10000×5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。 X=randn(10000,5); M=mean(X) D=std(X) R=corrcoef(X)

例6.6 对下列矩阵做各种排序。 185412613713-?? ??=?? ??-?? A A=[1,-8,5;4,12,6;13,7,-13]; sort(A) %对A 的每列按升序排序 -sort(-A,2) %对A 的每行按降序排序 [X,I]=sort(A) %对A 按列排序，并将每个元素所在行号送矩阵I 例6.7 给出概率积分 2 (d x x f x x -? e 的数据表如表6.1所示，用不同的插值方法计算f (0.472)。 x=0.46:0.01:0.49; %给出x ，f(x) f=[0.4846555,0.4937542,0.5027498,0.5116683]; format long interp1(x,f,0.472) %用默认方法，即线性插值方法计算f(x) interp1(x,f,0.472,'nearest') %用最近点插值方法计算f(x) interp1(x,f,0.472,'spline') %用3次样条插值方法计算f(x) interp1(x,f,0.472,'cubic') %用3次多项式插值方法计算f(x) format short 例6.8 某检测参数f 随时间t 的采样结果如表6.2，用数据插值法计算t =2，7，12，17，22，17，32，37，42，47，52，57时的f 值。 T=0:5:65; X=2:5:57;

数值计算方法与Matlab样卷答案

腹有诗书气自华 《数值计算方法与Matlab 》 样卷答案 一．填空题：（每空3分,共42分） 1． 8，6105.0-? 。 2．)(3)1(2)1(1)(3)1(2)1(1)(3)1(3)(3)(2)1(1)(3)(2)1(1)(2)1(2)(3)(2)(1)(3)(2)(1)(1)1(1)1(22)22()1()1(222)1()222(k k k k k k k k k k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x x x x x x x x x x ωωωωωωωωωω ωωωω-+--=---?+=+--+-=---?+=++--=+--?+=+++++++++， )2,1(∈ω。 3．],[1b a C S m -∈。4. 1e 2e ---x ，???==-=?--? ,3,2,1,0;0,e 1d )(e 110k k x x g k x ，正交投影。 5. 2阶，6阶。 6．10.6658，10.9521，10.9501。 7. 4002.2)00.1(=ε，4030.2)01.1(=ε。 二．解下列各题：（每题9分，共36分） 1．解：令)1(2 3+=t x ， （2分） 则??-+++=+1123 02 dt )1(25.21)1(49d 1t t x x x ???++++???++-+-≈22)6.01(25.21)6.01(9525.219 8)6.01(25.21)6.01(9549 （8分） 210631.10≈ （9分） 2．解：记系数矩阵为A, 对增广矩阵[]b A |作初等行运算, ??????????--401533933112??????????--==5.55.115 .35.405.75.401125.1,5.11,31,2l l ??????????---=45.114005.75.4011212,3l , 所以13-=x ,2)5.75.1(5.4112=-=x x ,1)1(2 1321=-+-=x x x ,即方程组的解为 [1,2,-1]T . （4分） 故系数矩阵A 的LU 分解为???? ??????--???????????---=4005.75.40112115.1015.1001A 。 （6分）

Matlab关于数值计算的实现

Matlab关于数值计算的实现 摘要：数值计算（numerical computation computation），主要研究更好的利用计算机更好的进行数值计算，解决各种数学问题。数值分析包括离散傅里叶变换，考虑截断误差，计算误差，函数的敛散性与稳定性等。在数学方面，数值计算的主要研究数值微分与积分，数据的处理与多项式计算，最优化问题，线性方程与非线性方程的求解，常微分方程的数值求解等。同时，数值计算在物理，化学，经济等方面也有研究，本文暂且不表。M atlab软件历经二十多年来的发展，已成为风靡世界的数学三大软件（matlb，Mathematica l，Maple）之一，在数学类科技应用软件中在数值计算方面首屈一指。Matlab以矩阵为数据操作的基本单位，使得矩阵运算十分便捷快速，同时Matlab还提供了海量的计算函数，而且使用可靠地算法进行计算，能使用户在繁复的数学运算中解脱，Matlab还具有方便且完善的图形处理功能，方便绘制二维和三维图形并修饰。

目录 1.数值分析(离散傅里叶变换，考虑截断误差，计算误差，函数 的敛散性与稳定性) 2.数值计算(数值微分与积分，数据的处理与多项式计算， 最优化问题，线性方程与非线性方程的求解，常微分方程的数值求解) 3.图形处理功能(方便绘制二维和三维图形并修饰) 4.总结

1.数据统计与分析 Matlab 可以进行求矩阵的最大最小元素，平均值与中值，关于矩阵元素的求和与求积，累加和与累乘积，标准方程，相关系数，元素排序。现在以求标准方差举例说明Matlab 的实现。 在Matlab 中，实现标准方差计算的函数为std 。对于向量(Y ），std （Y ）实现返回一个标准方差，而对于矩阵（A ），std （A ）返回一个行向量，该行向量的每个元素对应着矩阵A 各行或各列的标准方差。一般调用std 函数的格式为std （A ，flag ，dim ） Dim 取1或者2分别对应求各列或各行的标准方差，flag 取1时，按照标准方差的计算公式 ∑-=-=N i x x S i N 1 2 1)(11来计算。若flag 取2，则用公式 ∑-==N i x x S i N 1 2 2) (1 进行计算。默认的flag 取值为0，dim 取值为1。课本page143 2. 离散傅里叶变换 离散傅里叶变换广泛应用于信号的分析，光谱和声谱分析、全息技术等各个领域。但直接计算dft 的运算量与变化的长度N 的平方成正比，当N 较大时，计算量太大。随着计算机技术的迅速发展，在计算机上进行离散傅里叶变换计算成为可能。特别是快速傅里叶变换算法的出现，为傅里叶变换的应创造了条件。 （1）：傅里叶变换算法的简述。 傅立叶变换是一种分析信号的方法，它可分析信号的成分，也可用这些成分合成信号。许多波形可作为信号的成分，比如正弦波、方波、锯齿波等，傅立叶变换用正弦波作为信号的成分. f(t)是t 的周期函数，如果t 满足狄里赫莱条件：在一个以2T 为周期内f(X)连续或只有有限个第一类间断点，附f （x ）单调或可划分成有限个单调区间，则F （x ）以2T 为周期的傅里叶级数收敛，和函数S （x ）也是以2T 为周期的周期函数，且在这些间断点上，函数是有限值；在一个周期内具有有限个极值点；绝对可积。则有下图①式成立。称为积分运算f(t)的傅立叶变换， ②式的积分运算叫做F(ω)的傅立叶逆变换。F(ω)叫做f(t)的像函数，f(t)叫做 F(ω)的像原函数。F(ω)是f(t)的像。f(t)是F(ω)原像。 ①傅立叶变换 ②傅立叶逆变换

MATLAB实验二MATLAB的数值运算和程序

课程名称：Matlab语言 开设时间：2016—2017学年第 2 学期 专业班级：学生学号：学生姓名： 实验名称：实验二、MATLAB的数值运算和程序实验成绩： 指导教师：批改时间： 一、实验目的和要求 1）掌握基本的矩阵运算及常用的函数。 2）掌握MATLAB函数的编写及调试方法。 3）掌握MATLAB常用的数值运算函数。 二、实验仪器和设备 计算机一台 三、实验过程 1、一维数组在命令窗口执行下面指令，观察输出结果，体味数组创建和寻访方法，%号后面的为注释，不用输入。 rand('state',0) % 把均匀分布伪随机发生器置为0 状态 x=rand(1,5) % 产生（1*5）的均布随机数组 x(3) % 寻访数组x 的第三个元素。 x([1 2 5]) % 寻访数组x 的第一、二、五个元素组成的子数组。 x(1:3) % 寻访前三个元素组成的子数组 x(3:end) % 寻访除前2 个元素外的全部其他元素。end 是最后一个元素的下标。 x(3:-1:1) % 由前三个元素倒排构成的子数组 x(find(x>0.5)) % 由大于0.5 的元素构成的子数组 x([1 2 3 4 4 3 2 1]) % 对元素可以重复寻访，使所得数组长度允许大于原数组。 x(3) = 0 % 把上例中的第三个元素重新赋值为0。 x([1 4])=[1 1] % 把当前x 数组的第一、四个元素都赋值为1。 x[3]=[] % 空数组的赋值操作

2、在命令窗口执行下面指令，观察输出结果 a=2.7358; b=33/79; % 这两条指令分别给变量 a ， b 赋值。 C=[1,2*a+i*b,b*sqrt(a);sin(pi/4),a+5*b,3.5+i] % 这指令用于创建二维组C M_r=[1,2,3;4,5,6],M_i=[11,12,13;14,15,16] % 创建复数数组的另一种方法 CN=M_r+i*M_i % 由实部、虚部数组构成复数数组 3. 记录下面题目的程序和运行后的结果。 1??????=654321a ??????-=531142b ???? ? ?????-=201c ??????????=063258741d 下列运算是否合法，为什么？如合法，结果是多少？

数值分析 matlab 实验4

(1) 解题过程如下： （1）MATLAB中创建复化梯形公式和复化辛普森公式的 M 文件：1）复化梯形公式文件： function s=T_fuhua(f,a,b,n) h=(b-a)/n; s=0; for k=1:(n-1) x=a+h*k; s=s+feval(f,x); end s=h*(feval(f,a)+feval(f,b))/2+h*s; 2）复化辛普森公式文件： function s=S_fuhua(f,a,b,n) h=0; h=(b-a)./(2*n); s1=0; https://www.wendangku.net/doc/fd16545220.html, -5- s2=0; for k=1:n-1 x=a+h*2*k; s1=s1+feval(f,x); end for k=1:n x=a+h*(2*k-1); s2=s2+feval(f,x); end

s=h*(feval(f,a)+feval(f,b)+s1*2+s2*4)/3; 在MATLAB中输入： f=inline('x/(4+x^2)');a=0;b=1; %inline 构造内联函数对象 for n=2:10 s(n-1)=T_fuhua(f,a,b,n);s(n-1)=vpa(s(n-1),10); %调用复化梯形公式，生成任意精度的数值 end exact=int('x/(4+x^2)',0,1);exact=vpa(exact,10) %求出积分的精确值 输出结果：exact = .1115717755 s = Columns 1 through 6 0.1088 0.1104 0.1109 0.1111 0.1113 0.1114 Columns 7 through 9 0.1114 0.1114 0.1115 在MATLAB中输入以下函数用以画出计算误差与 n 之间的曲线： r=abs(exact-s); n=2:10; plot(double(n),double(r(n-1))) 得到结果如图所示： （2）在 MATLAB中输入以下程序代码： f=inline('x/(4+x^2)');a=0;b=1;n=9; %inline 构造内联函数对象 t=T_fuhua(f,a,b,n);t=vpa(t,10) s=S_fuhua(f,a,b,n);s=vpa(s,10)

数值分析实验— MATLAB实现

数值分析实验 ——MATLAB实现 姓名sumnat 学号2013326600000 班级13级应用数学2班 指导老师 2016年1月

一、插值：拉格朗日插值 (1) 1、代码： (1) 2、示例： (1) 二、函数逼近：最佳平方逼近 (2) 1、代码： (2) 2、示例： (2) 三、数值积分：非反常积分的Romberg算法 (3) 1、代码： (3) 2、示例： (4) 四、数值微分：5点法 (5) 1、代码： (5) 2、示例： (6) 五、常微分方程：四阶龙格库塔及Adams加速法 (6) 1、代码：四阶龙格库塔 (6) 2、示例： (7) 3、代码：Adams加速法 (7) 4、示例： (8) 六、方程求根：Aitken 迭代 (8) 1、代码： (8) 2、示例： (9) 七、线性方程组直接法：三角分解 (9) 1、代码： (9) 2、示例： (10) 八、线性方程组迭代法：Jacobi法及G-S法 (11) 1、代码：Jacobi法 (11) 2、示例： (12) 3、代码：G-S法 (12) 4、示例： (12) 九、矩阵的特征值及特征向量：幂法 (13) 1、代码： (13) 2、示例： (13)

一、插值：拉格朗日插值 1、代码： function z=LGIP(x,y)%拉格朗日插值 n=size(x); n=n(2);%计算点的个数 syms a; u=0;%拉格朗日多项式 f=1;%插值基函数 for i=1:n for j=1:n if j==i f=f; else f=f*(a-x(j))/(x(i)-x(j)); end end u=u+y(i)*f;f=1; end z=expand(u);%展开 2、示例： >> x=1:6; y1=x.^5+3*x.^2-6; y2=sin(x)+sqrt(x); >> f1=LGIP(x,y1) f1 = -6+3*a^2+a^5 %可知多项式吻合得很好 >> f2=vpa(LGIP(x,y2),3) f2 = .962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5

数值计算方法实验指导(Matlab版)

数值计算方法实验指导(Matlab 版)

《数值计算方法》 实验指导 （Matlab版） 肇庆学院数学与统计学学院 计算方法课程组

《数值计算方法》实验1报告 班级： 20xx 级XXXXx 班 学 号： 20xx2409xxxx 姓名： XXX 成绩： 1. 实验名称 实验1 算法设计原则验证（之相近数相减、大数吃小数和简化计算步骤） 2. 实验题目 (1) 取16 10=z ，计算 z z -+1和) 1/( 1z z ++，验证 两个相近的数相减会造成有效数字的损失． (2) 按不同顺序求一个较大的数（123）与1000个较小的数（15 310-?）的和，验证大数吃小 数的现象． (3) 分别用直接法和秦九韶算法计算多项式 n n n n a x a x a x a x P ++++=--1 1 1 )(Λ 在x =1.00037处的值．验证简化计算步骤能减少运算时间． 对于第(3)题中的多项式P (x )，直接逐项计算 需要21 12)1(+=+++-+n n n Λ次乘法和n 次加法，使用秦

九韶算法 n n a x a x a x a x a x P ++++=-)))((()(1210Λ 则只需要n 次乘法和n 次加法． 3. 实验目的 验证数值算法需遵循的若干规则． 4. 基础理论 设计数值算法时，应避免两个相近的数相减、防止大数吃小数、简化计算步骤减少运算次数以减少运算时间并降低舍入误差的积累．两相近的数相减会损失有效数字的个数，用一个大数依次加小数，小数会被大数吃掉，乘法运算次数太多会增加运算时间． 5. 实验环境 操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程 (1) 直接计算并比较； (2) 法1：大数逐个加1000个小数，法2：先把1000个小数相加再与大数加； (3) 将由高次项到低次项的系数保存到数组A[n]中，其中n 为多项式次数．

数值分析五个题目的C语言及Matlab程序

本文档包含上一个文档中的五个数值分析实验题C语言程序及Matlab程序实验一 C程序 #include "stdio.h" #include "math.h" void main() { inti=0; float a=0.1,b=1.9,t=0.0,e=1.9; if((pow(a,7)-28*pow(a,4)+14)*(pow(b,7)-28*pow(b,4)+14)<0) if((7*pow(x,6)-112*pow(x,3))) printf("x=%f,i=%d,e=%f\n",x,i,e); for(i=1;i<7&&e>0.00001;i++) { t=x; x=x-(pow(x,7)-28*pow(x,4)+14)/(7*pow(x,6)-112*pow(x,3)); e=fabs(t-x); printf("x=%f,i=%d,e=%f\n",x,i,e);

} } Matable 程序 i=0; x=1.9;t=0.0;e=1.9; disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]); for i=1:7 t=x; x=x-(x^7-28*x^4+14)/(7*x^6-112*x^3); e=abs(t-x); disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]); if e<0.00001 break; end end 实验二 C程序 #include"stdio.h" #include"math.h" //已知量 double x[10]={1,2,3,4,5,6,7,8,9,10};

数值计算方法matlab程序

function [x0,k]=bisect1(fun1,a,b,ep) if nargin<4 ep=1e-5; end fa=feval(fun1,a); fb=feval(fun1,b); if fa*fb>0 x0=[fa,fb]; k=0; return; end k=1; while abs(b-a)/2>ep x=(a+b)/2; fx=feval(fun1,x); if fx*fa<0 b=x; fb=fx; else a=x; fa=fx; k=k+1; end end x0=(a+b)/2; >> fun1=inline('x^3-x-1'); >> [x0,k]=bisect1(fun1,1.3,1.4,1e-4) x0 = 1.3247 k = 7 >> 简单迭代法 function [x0,k]=iterate1(fun1,x0,ep,N) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep;

while abs(x-x0)>ep & k> fun1=inline('(x+1)^(1/3)'); >> [x0,k]=iterate1(fun1,1.5) x0 = 1.3247 k = 7 >> fun1=inline('x^3-1'); >> [x0,k]=iterate1(fun1,1.5) x0 = Inf k = 9 >> Steffesen加速迭代（简单迭代法的加速）function [x0,k]=steffesen1(fun1,x0,ep,N) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep; k=0; while abs(x-x0)>ep & k

解非线性方程的数值计算方法 用Matlab实现

湖北民族学院《数值分析》实验报告 实验名称 解非线性方程的数值计算方法 实验时间 2011年 11月 9日 姓名 王亚敏 班级 0209408 学号 020940807 成绩 实验报告内容要求： 一、实验目的与要求；二、实验内容；三、算法描述(数学原理或设计思路、计算公式、计算步骤)； 四、程序代码；五、数值结果；六、计算结果分析（如初值对结果的影响；不同方法的比较；该方法的特点和改进等）；七、实验中出现的问题，解决方法及体会（整个实验过程中（包括程序编写，上机调试等）出现的问题及其处理等广泛的问题）. 一、实验目的： ① 掌握二分法、不动点迭代法、牛顿迭代法、Steffensen 迭代法等常用的非线性方程迭代算法； ② 培养编程与上机调试能力. 二、实验内容： 1．请分别用二分法、不动点迭代法、牛顿迭代法、Steffensen 迭代法求方程 32 ()330f x x x x =+--=在1.5 附近的根. 参考答案：原方程在1.5 附近的根为1.732051x = 2．请自选两个非线性方程，分别用二分法、不动点迭代法、牛顿迭代法、Steffensen 迭代法求解. 三、算法 （一） 1．二分法 计算()0f x =的二分法如下： ① 输入求根取间[,]a b 和误差控制量ε，定义函数. 如果 ()()0f a f b <，转②；否则退出选用其它求根方法 ② 当a b ε->时，计算中点()/2x a b =+以及的值； 分情况处理 ()f x ε<：停止计算，，转④ ()()0f a f x <：修正区间[,][,]ax ab → ()()0f x f b <：修正区间[,][,]xb ab → ③ *2a b x += ④ 输出近似根* x 2．不动点(简单)迭代法算法思想 先将f(x)=0化成x=φ(x),然后取x 0,计算 x n+1=φ(x n ) n=0,1,2,… 若要求x * 满足f(x * )=0,则x * =)(*x ?，称x * 为函数)(x ?的一个不动点。