数值分析报告

重庆交通大学计算机与信息学院

设计性实

验报告

班 级： 05计科(3)班

实验项目性质： 综合设计性实验

实验所属课程： 数 值 分 析

实验室(中心) ： 计算机软件中心

指 导 教 师 ： 程 攀

学 生 学 号 ： 05060315

学 生 姓 名 ： 郭 盛

实验完成时间： 2008 年 1 月 14 日

实验一 非线性方程求根

一. 实验目的:

1. 学会用Newton 迭代法求非线性方程的根

2. 通过实际编程计算进一步了解该方法的优缺点

二. 实验学时:

课内实验:2课时

三. 实验算法:

实验算法(Newton算法):

（1） 取初始点x 0，最大迭代次数N 和精度要求ε，置k=0。

（2） 如果f ’(X k )=0,则停止计算，否则计算

Xk+1=X k -f(X k )/f’(X k )

（3） 若| Xk+1- X k |

（4） 若k=N，则停止计算；否则置k=k+1，转（2）。

四. 实验内容:

用迭代法求方程e +10*x-2=0的近似根，要求误差不超过0.5*10，取初值x 0=0。 -3x

五. 实验代码及结果:

newton 迭代：

function [p1,err,k,y]=newton(f,df,p0,delta,N)

%f是非线性函数

%df是f 的微商

%p0是初始值

%delta是给定允许误差

%N迭代的最大次数

%p1是牛顿法求得的方程的近似值

%err是p0的误差估计

%k是迭代次数

%y=f(p1)

p0,feval('f',p0)

for k=1:N

p1=p0-feval('f',p0)/feval('df',p0);

err=abs(p1-p0);

p0=p1;

p1,err,k,y=feval('f',p1)

if(err

break,end

p1,err,k,y=feval('f',p1)

end

方程函数：

function y=f(x)

y=exp(x)+10*x-2

微分函数：

function y=df(x)

y=exp(x)+10;

实验结果：

newton('f','df',0,5^(-4),10)

p0 =0

y =-1

ans =-1

y =-1

p1 =0.0909

err =0.0909

k =1

y =0.0043

y =0.0043

p1 =0.0909

err =0.0909

k =1

y =0.0043

y =0.0043

y =0.0043

p1 =0.0905

err =3.8398e-004

k =2

y =8.0727e-008

y = 8.0727e-008

ans =0.0905

六. 实验总结:

通过本次实验, 加深了对非线性问题的迭代法与线性问题迭代法的差别的认识，深入了解了迭代法及初始值与迭代收敛性的关系。迭代法是求解非线性方程的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。

实验二 线性代数方程组的解法

--------列主元Gauss 消元法

一. 实验目的

1. 学会用列主元Gauss 消去法求解线性代数方程组

2. 加深对选取主元素的重要性的认识

3. 通过实际变成计算进一步了解该方法的优缺点

二. 实验算法:

实验算法（列主元Gauss 消去法）

（1） 输入数据A 和b

（2） 对于k=1,2,…..,n-1, 做

A, 按列选主元。令|a r k|=max |a i k|

K

如果| a r k |=0，则停止计算

B, 交换两行。如果r>k，则a k ja r j j=k,k+1,…..,n b k b r

C, 消元计算。置m i k= a i k /a k k i=k+1,k+2,…,n

a i j=a i j-m i k*a k j i,j=k+1,k+2,…,n

b i = bi -m i k* b k i=k+1,k+2,…n

(3) 如果a n n=0，则停止计算，否则进行回代，对于k=n,n-1,…1，置

x k =(b k -∑a k j*xj ) /akk （当下标大于上标时，不做求和运算）

(4)输出x 1,x 2, …,x n ，即为方程组Ax=b的解

三. 实验内容:

用列主元Gauss 消去法求解线性代数方程组

1.1348 3.8326 1.1651 3.4017

x1 9.5342

0.5301 1.7875 2.5330 1.5435 x2 6.3941

3.4129 4.9317 8.7643 1.3142 x3 = 1.2371 4.9998 10.6721 0.0147

四. 实验代码及结果:

function u = gauss (a, b)

n = length (b) ;

%success = 1; % 判断gauss 消去法是否成功，1成功，0失败

b_tmp = 1:n; % 记录列变换的情况

% 消去过程begin

for k=1: n-1

% 选主元begin

max = abs(a(k, k)); % 当前主元的最大值

maxL = k; % 当前主元所在的行

maxC = k; % 当前主元所在的列

for j = k:n

for i = k:n

if abs(a(j, i)) > max

%[maxL, maxC] = [j, i];

maxL = j;

maxC = i;

max = abs(a(j, i));

end

end

end

% 交换第k 行和第maxL 行begin

for j = k:n

tmp = a(maxL, j);

a(maxL, j) = a(k, j);

a(k, j) = tmp;

end

tmp = b(maxL);

b(maxL) = b(k);

b(k) = tmp;

% 交换第k 行和第maxL 行end

% 交换第k 列和第maxC 列begin

for j = 1:n

tmp = a(j, maxC);

a(j, maxC) = a(j, k);

a(j, k) = tmp;

end

tmp = b_tmp(maxC);

b_tmp(maxC) = b_tmp(k);

b_tmp(k) = tmp;

% 交换第k 列和第maxC 列 end

% 选主元end

for i = k+1 : n

if abs (a( k, k)) > 1e-6

mult = a(i, k) / a(k, k) ;

for j =k+1 : n

a(i, j) = a(i, j) - mult * a(k, j) ;

end % end for

b(i) = b(i) - mult * b(k) ;

else

%success = 0; % Gauss消去法失败

disp ( '列主元Gauss 消去法失败');

pause;

exit;

end % end if

end % end for

end % end for

% 消去过程end

% 回代过程begin

x(n)= b(n) / a(n, n) ;

for i = n-1:-1:1

s=0 ;

for j = i+1: n

s=s+a(i, j) * x(j) ;

end

x(i) = (b(i) - s) / a(i, i) ;

end

% 回代过程end

% 获得正确的解

for i = 1:n

u(b_tmp(i)) = x(i);

end

End

五. 实验总结:

通过本次实验基本上掌握并学会用列主元Gauss 消去法求解线性代数方程组同时也加深了对选取主元素的重要性的认识。Gauss 消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何确保Gauss 消去法作为数值算法的稳定性呢才是关键,Gauss 消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验三 线性代数方程组的解法

——Gauss-Seidel 迭代法

一. 实验目的:

1. 学会用Gauss-Seidel 迭代法求解线性代数方程组；

2. 进一步领会Gauss-Seidel 法的迭代思想；

3. 通过实际编程计算进一步了解该方法的优缺点。

二. 实验算法:

实验算法：（Gauss-Seidel 迭代法）：

（1） 输入数据A 和b 。

（2） 取初始点x (0)，置K=0，精度要求和最大迭代次数N 。

（3） 按下式计算x (k+1),

X i (k+1)=(1/aii )(bi -∑a ij x j (k+1)-∑a ij x j (k)) i=1,2,…,n.

(4) 若||x(k+1)-x (k)||2

（5） 若K=N，则停止计算（输出某些信息）；否则置K=K+1，转（3）。

三. 实验内容:

例 用Gauss-Seidel ||x(k+1)-x (k)||2

-1 4 -1 0 -1 0 x2 5

0 -1 4 -1 0 -1 x3 = -2

-1 0 -1 4 -1 0 x4 5

0 -1 0 -1 4 -1 x5 -2

0 0 -1 0 -1 4 x6 6

四. 实验代码：

function[x,k,flag]=Gau_Seid(A,b,delta,max1)

%Gauss-Seidel迭代法

%A为方程组的系数矩阵

%b为方程组的右端

%delat为精度要求，实验中为0.0001

%max1为最大迭代次数，实验中设为50

%x为方程组的解;%k为迭代次数;

%flag为指标变量 flag='ok！' 表示迭代收敛到指标要求 ；flag='fail'表示迭代失败 if nargin

if nargin

n=length(A);k=0;

x=zeros(n,1);y=zeros(n,1);flag='ok';

while 1

y=x;

for i=1;n

z=b(i);

for j=1:n

if j~=i

z=z-A(i,j)*x(j);

end

end

if (abs(A(i,j))

flag='fail';return;

end

z=z/A(i,i);x(i)=z;

end

if norm(y-x,inf)

break;

end

k=k+1;

end

五. 实验结果:

>> A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0 ;0 -1 4 -1 0 -1 ;

-1 0 -1 4 -1 0;0 -1 0 -1 4 -1; 0 0 -1 0 -1 4]

A =

4 -1 0 -1 0 0

-1 4 -1 0 -1 0

0 -1 4 -1 0 -1

-1 0 -1 4 -1 0

0 -1 0 -1 4 -1

0 0 -1 0 -1 4

>> b=[0 5 -2 5 -2 6]

b =0 5 -2 5 -2 6

>> delta=0.0001

delta =1.0000e-004

K>> Gau_Seid(A,b,delta,max1)

11 if nargin

K>> Gau_Seid(A,b,delta,max1)

11 if nargin

六. 实验总结:

通过这次实验, 可分析如下结果:

1) 优点分析:

1. 程序容错性好，对于能够求解和不能够求解的方程程序都会进行处理并给予回复.

2. 程序简单易懂，可读性强。

3. 执行效率高，即使是大型矩阵也可以很快求解

2)缺点分析:

1. 程序每次只能对一个方程组进行求解，对于矩阵A ，b ，程序能一次调用就全部求解。

2. 算法没能很好的利用MATLAB 强大的矩阵运算功能来减少代码，使得代码冗长。

实验四 插值法

————三次样条插值法

一. 实验目的:

1． 学会用三次样条函数作插值计算

2． 通过应用和实际编程计算进一步了解该方法的原理和优缺点。

二. 实验算法:

实验算法（求三次样条插值勤函数S （x ）在一列点x 1,x 2, „x m 上的近似值）；

（1） 输入x i ,y i ,(i=0,1,„,n), 边界条件y 0,y n , 及x j (j=1,2,„m).

（2） 计算h i .h i =xi ,x i-1 i=1,2,„n.

（3） 按下列公式计算μi , λi,g i , i=1,2,„n-1

μi =hi /(hi +hi+1), λi , μi

g i =6/(hi +hi+1)(yi+1-y i )/hi+1-(yi -y i-1)/hi )

(4) 按下列公式计算g 0,g n ,

g 0=(6/h1)((y1-y 0)/h1-y 0`)

g n =6/hn (yn `-(yn -y n-1/hn ))

(5) 用追赶法求解如下方程组得M i , i=0,1,„n,

0 0

μ1 2 λ1 11

„ „ „﹕

μn-1 2 λg n-1

n n

（6）确定点x j (j=1,2,„m) 所在的区间[xi-1,xi],并按下列公式计算S （xj ）; S （x j ）3322=Mi-1(xi -x) +Mi ((x-xi-1) /6hi )+(yi-1-(Mi-1/6)hi )((xi -x)/hi )+(yi -(Mi /6)hi )((x-xi-1)/hi ) (x∈[xi-1,xi],i=1,2,„n)

三. 实验内容:

已知直升飞机旋转机翼外开曲线轮廓线（如图）上的某些弄值点（见表）及端点处的一阶导数数值

y`(x0)=1.86548, y(x18)=-0.046115

试计算该曲线上横坐标为

2，4，6，12，16，30，60，110，180，280，280，400，515

四. 实验代码:

function S=csfit(x,y,dx0,dxn)

%x,y分别为n 个接点的横坐标所组成的向量及纵坐标组成的向量

%dx(),dn分别为S 的导数在x0,xn 处的值

n=length(x)-1;

h=diff(x);

d=diff(y)./h;

a=h(2:n-1);

b=2*(h(1:n-1)+h(2:n));

c=h(2:n);

u=6*diff(d);

b(1)=b(1)-h(1)/2;

u(1)=u(1)-3*(d(1)-dx0);

for k=2:n-1

temp=a(k-1)/b(k-1);

b(k)=b(k)-temp*c(k-1);

u(k)=u(k)-temp*u(k-1);

end

m(n)=u(n-1)/b(n-1);

for k=n-2:-1:1

m(k+1)=(u(k)-c(k)*m(k+2))/b(k);

end

m(1)=3*(d(1)-dx0)/h(1)-m(2)/2;

m(n+1)=3*(dxn-d(n))/h(n)-m(n)/2;

for k=0:n-1

S(k+1,1)=(m(k+2)-m(k+1))/(6*h(k+1));

S(k+1,2)=m(k+1)/2;

S(k+1,3)=d(k+1)-h(k+1)*(2*m(k+1)+m(k+2))/6;

S(k+1,4)=y(k+1);

End

五. 实验结果

>> x=[0,1,2,3];y=[0,0.5,2,1.5];dx0=0.2;dxn=-1;

>> s=csfit(x,y,dx0,dxn)

s =

0.4731 -0.1731 0.2000 0

-1.0192 1.2462 1.2731 0.5000

0.6558 -1.8115 0.6558 2.0000

六. 实验总结:

通过本次实验, 考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，L n (x ) 是否也更加靠近被逼近的函数。龙格（Runge ）给出一个例子是极著名并富有启发性的。项式插值是不收敛的，即插值的节点多，效果不一定就好。

按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数, 可以用MATLAB 的函数“spline ”作此函数的三次样条插值。

实验五 数值积分法

一. 实验目的:

1. 学会用复合辛普森公式计算某积分的近似值；

2. 学会用半步长复合梯形法计算某积分的近似值；

3. 通过实际编程计算进一步了解方法的原理和优缺点。

二. 实验算法:

实验算法：

1. 复合辛普森公式

∫a f（x ）≈h/6[f(a)+4∑f(Xk+1/2)+2∑f(Xk )+f(b)]

其中X k =a+k*h,(k=0，1，...,n),h=（b-a ）/n,Xk+1/2=Xk +h/2。

2. 半步长复合梯形法的递推算法：

T1=（b-a ）/2*[f(a)+f(b)]

T2=1/2*Tn+h/2∑f(Xk+1/2)

其中X k =a+k*h, (k=0，1，...,n-1), h=（b-a ）/n,Xk+1/2=Xk +h/2。 b

三. 实验内容：

1. 用复合辛普森公式计算积分∫0 48（1+cos2x ）1/2 dx ，是误差不超过10-4

2. 用半步长复合梯形法的递推算法计算积分

∫01（1-e -x ）/x dx

使误差不超过10-4

四. 实验代码:

function s=trapr1(f,a,b,n)

h=(b-a)/n;

s=0;

for k=1:n-1

x=a+h*k;

s=s+feval('lp',x);

end

s=h*(feval('lp',a)+feval('lp',b))/2+h*s;

function z=lp(x)

z=1/(1+x^2);

function z=lp(x)

z=(1-exp(-x)^(1/2))/x;

五. 实验结果:

>> trapr1('lp',0.00000001,1,10)

ans =

0.4439

>>

>> trapr1('lp',-1,1,10)

ans =

1.5675

六. 实验总结:

通过本次实验可知, 线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。

高维空间中的积分，如果维数不很高且积分区域是规则的或者能等价地写成多重积分的形式，可以用一元函数积分的数值方法来计算高维空间的积分。蒙特卡罗方法对计算复杂区域甚至不连通的区域上的积分并没有特殊的困难。

实验六 常微分方程的数值解法

一. 实验目的:

1. 学会用改进的欧拉法求一阶常微分方程初值问题的数值解；

2. 学会用经典R-K 法求一阶常微分方程初值问题的数值解；

3. 通过编程计算进一步了解方法的优缺点

二. 实验算法:

实验算法：

1. 改进的欧拉公式：

T 1=y i +h*f(x i ,y i )

T 2=y i +h*f(x i +1,T1)

y i+1=(T1+T2)/2

2. 经典R-K 公式：

y i+1=y i +h/6*(K1+2*K2+2*K3+K4)

K 1=f(x i ,y i )

K 2=f(x i +h/2,y i +h/2*K1)

K 3=f(x i +h/2,y i +h/2*K2)

K 4=f(x i +h,y i+h*K3)

三. 实验内容：

用下列方法求

-2 dy/dx=2/3*x*y,x ∈[0,1]

y(0)=1

21/3的数值解（取h=0.1），并将计算结果与准确解y=(1+x) 进行比较：

（1） 改进的欧拉法

（2） 经典R-K 法

四. 实验代码:

function E=Euler(fun,x0,y0,xN,N)

x=zeros(1,N+1);y=zeros(1,N+1);

x(1)=x0;y(1)=y0;

h=(xN-x0)/N;

for n=1:N

x(n+1)=x(n)+h;

z0=y(n)+h*feval(fun,x(n),y(n));

y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0));

end

T=[x',y']

function z=f(x,y)

z=x+y;

function R=rk4(f,a,b,ya,N)

h=(b-a)/N;

T=zeros(1,N+1);

Y=zeros(1,N+1);

T=a:h:b;

Y(1)=ya;

for j=1:N

k1=h*feval(f,T(j),Y(j));

k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);

k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);

k4=h*feval(f,T(j)+h,Y(j)+k3);

Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;

end

R=[T' Y']

function rk=rk(x,y)

rk=y-2*x/y;

五. 实验结果:

%取迭代次数为5

>> Euler('f',0,1,0.1,5)

T =

0 1.0000

0.0200 1.0204

0.0400 1.0416

0.0600 1.0637

0.0800 1.0866

0.1000 1.1103

function a=b(x,y)

a=(2/3)*x/y^2;

>> Euler('b',0,1,0.1,5)

T =

0 1.0000

0.0200 1.0001

0.0400 1.0005

0.0600 1.0012

0.0800 1.0021

0.1000 1.0033

%R-K法

%迭代次数为10

>> rk4('rk',0,1,1,10)

R =

0 1.0000

0.1000 1.0954

0.2000 1.1832

0.3000 1.2649

0.4000 1.3416

0.5000 1.4142

0.6000 1.4832

0.7000 1.5492

0.8000 1.6125

0.9000 1.6733

1.0000 1.7321

ans =

0 1.0000

0.1000 1.0954

0.2000 1.1832

0.3000 1.2649

0.4000 1.3416

0.5000 1.4142

0.6000 1.4832

0.7000 1.5492

0.8000 1.6125

0.9000 1.6733

1.0000 1.7321

实验

>> rk4('b',0,1,1,10)

R =

0 1.0000

0.1000 1.0033

0.2000 1.0132

0.3000 1.0291

0.4000 1.0507

0.5000 1.0772

0.6000 1.1079

0.7000 1.1422

0.8000 1.1793

0.9000 1.2187

1.0000 1.2599

ans =

0 1.0000

0.1000 1.0033

0.2000 1.0132

0.3000 1.0291

0.4000 1.0507

0.5000 1.0772

0.6000 1.1079

0.7000 1.1422

0.8000 1.1793

0.9000 1.2187

1.0000 1.2599

六. 实验总结:

通过本次实验, 可得到如下结果：

[结果分析]：

从下面的图来看，我们虽然用的是不同的方法，但是求出来的结果却是非常吻合。

[程序的优缺点]：

优点：函数通用性好（特别是经典4阶的Runge-Kutta 方法），只需要稍作修改就可以直接用到其他的微分方程组上。

缺点：简单Euler 方程精度不是很高。

改进的欧拉法也称为梯形法。不难发现，欧拉公式是关于y n +1的显式，即只要已知y n , 经过一次计算便可得y n +1的值，而改进的 欧拉公式是以y n +1的隐式方程给出，不能直接得到y n +1。隐式方程通常用迭代法求解，而迭代过程的实质是逐步显式化。

R-K 方法不是通过求导数的方法构造近似公式，而是通过计算 不同点上的函数值，并对这些函数值作线性组合，构造近似公式，再把近似公式与解的泰勒展开 式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。

重庆交通大学计算机与信息学院

设计性实

验报告

班 级： 05计科(3)班

实验项目性质： 综合设计性实验

实验所属课程： 数 值 分 析

实验室(中心) ： 计算机软件中心

指 导 教 师 ： 程 攀

学 生 学 号 ： 05060315

学 生 姓 名 ： 郭 盛

实验完成时间： 2008 年 1 月 14 日

实验一 非线性方程求根

一. 实验目的:

1. 学会用Newton 迭代法求非线性方程的根

2. 通过实际编程计算进一步了解该方法的优缺点

二. 实验学时:

课内实验:2课时

三. 实验算法:

实验算法(Newton算法):

（1） 取初始点x 0，最大迭代次数N 和精度要求ε，置k=0。

（2） 如果f ’(X k )=0,则停止计算，否则计算

Xk+1=X k -f(X k )/f’(X k )

（3） 若| Xk+1- X k |

（4） 若k=N，则停止计算；否则置k=k+1，转（2）。

四. 实验内容:

用迭代法求方程e +10*x-2=0的近似根，要求误差不超过0.5*10，取初值x 0=0。 -3x

五. 实验代码及结果:

newton 迭代：

function [p1,err,k,y]=newton(f,df,p0,delta,N)

%f是非线性函数

%df是f 的微商

%p0是初始值

%delta是给定允许误差

%N迭代的最大次数

%p1是牛顿法求得的方程的近似值

%err是p0的误差估计

%k是迭代次数

%y=f(p1)

p0,feval('f',p0)

for k=1:N

p1=p0-feval('f',p0)/feval('df',p0);

err=abs(p1-p0);

p0=p1;

p1,err,k,y=feval('f',p1)

if(err

break,end

p1,err,k,y=feval('f',p1)

end

方程函数：

function y=f(x)

y=exp(x)+10*x-2

微分函数：

function y=df(x)

y=exp(x)+10;

实验结果：

newton('f','df',0,5^(-4),10)

p0 =0

y =-1

ans =-1

y =-1

p1 =0.0909

err =0.0909

k =1

y =0.0043

y =0.0043

p1 =0.0909

err =0.0909

k =1

y =0.0043

y =0.0043

y =0.0043

p1 =0.0905

err =3.8398e-004

k =2

y =8.0727e-008

y = 8.0727e-008

ans =0.0905

六. 实验总结:

通过本次实验, 加深了对非线性问题的迭代法与线性问题迭代法的差别的认识，深入了解了迭代法及初始值与迭代收敛性的关系。迭代法是求解非线性方程的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。

实验二 线性代数方程组的解法

--------列主元Gauss 消元法

一. 实验目的

1. 学会用列主元Gauss 消去法求解线性代数方程组

2. 加深对选取主元素的重要性的认识

3. 通过实际变成计算进一步了解该方法的优缺点

二. 实验算法:

实验算法（列主元Gauss 消去法）

（1） 输入数据A 和b

（2） 对于k=1,2,…..,n-1, 做

A, 按列选主元。令|a r k|=max |a i k|

K

如果| a r k |=0，则停止计算

B, 交换两行。如果r>k，则a k ja r j j=k,k+1,…..,n b k b r

C, 消元计算。置m i k= a i k /a k k i=k+1,k+2,…,n

a i j=a i j-m i k*a k j i,j=k+1,k+2,…,n

b i = bi -m i k* b k i=k+1,k+2,…n

(3) 如果a n n=0，则停止计算，否则进行回代，对于k=n,n-1,…1，置

x k =(b k -∑a k j*xj ) /akk （当下标大于上标时，不做求和运算）

(4)输出x 1,x 2, …,x n ，即为方程组Ax=b的解

三. 实验内容:

用列主元Gauss 消去法求解线性代数方程组

1.1348 3.8326 1.1651 3.4017

x1 9.5342

0.5301 1.7875 2.5330 1.5435 x2 6.3941

3.4129 4.9317 8.7643 1.3142 x3 = 1.2371 4.9998 10.6721 0.0147

四. 实验代码及结果:

function u = gauss (a, b)

n = length (b) ;

%success = 1; % 判断gauss 消去法是否成功，1成功，0失败

b_tmp = 1:n; % 记录列变换的情况

% 消去过程begin

for k=1: n-1

% 选主元begin

max = abs(a(k, k)); % 当前主元的最大值

maxL = k; % 当前主元所在的行

maxC = k; % 当前主元所在的列

for j = k:n

for i = k:n

if abs(a(j, i)) > max

%[maxL, maxC] = [j, i];

maxL = j;

maxC = i;

max = abs(a(j, i));

end

end

end

% 交换第k 行和第maxL 行begin

for j = k:n

tmp = a(maxL, j);

a(maxL, j) = a(k, j);

a(k, j) = tmp;

end

tmp = b(maxL);

b(maxL) = b(k);

b(k) = tmp;

% 交换第k 行和第maxL 行end

% 交换第k 列和第maxC 列begin

for j = 1:n

tmp = a(j, maxC);

a(j, maxC) = a(j, k);

a(j, k) = tmp;

end

tmp = b_tmp(maxC);

b_tmp(maxC) = b_tmp(k);

b_tmp(k) = tmp;

% 交换第k 列和第maxC 列 end

% 选主元end

for i = k+1 : n

if abs (a( k, k)) > 1e-6

mult = a(i, k) / a(k, k) ;

for j =k+1 : n

a(i, j) = a(i, j) - mult * a(k, j) ;

end % end for

b(i) = b(i) - mult * b(k) ;

else

%success = 0; % Gauss消去法失败

disp ( '列主元Gauss 消去法失败');

pause;

exit;

end % end if

end % end for

end % end for

% 消去过程end

% 回代过程begin

x(n)= b(n) / a(n, n) ;

for i = n-1:-1:1

s=0 ;

for j = i+1: n

s=s+a(i, j) * x(j) ;

end

x(i) = (b(i) - s) / a(i, i) ;

end

% 回代过程end

% 获得正确的解

for i = 1:n

u(b_tmp(i)) = x(i);

end

End

五. 实验总结:

通过本次实验基本上掌握并学会用列主元Gauss 消去法求解线性代数方程组同时也加深了对选取主元素的重要性的认识。Gauss 消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何确保Gauss 消去法作为数值算法的稳定性呢才是关键,Gauss 消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验三 线性代数方程组的解法

——Gauss-Seidel 迭代法

一. 实验目的:

1. 学会用Gauss-Seidel 迭代法求解线性代数方程组；

2. 进一步领会Gauss-Seidel 法的迭代思想；

3. 通过实际编程计算进一步了解该方法的优缺点。

二. 实验算法:

实验算法：（Gauss-Seidel 迭代法）：

（1） 输入数据A 和b 。

（2） 取初始点x (0)，置K=0，精度要求和最大迭代次数N 。

（3） 按下式计算x (k+1),

X i (k+1)=(1/aii )(bi -∑a ij x j (k+1)-∑a ij x j (k)) i=1,2,…,n.

(4) 若||x(k+1)-x (k)||2

（5） 若K=N，则停止计算（输出某些信息）；否则置K=K+1，转（3）。

三. 实验内容:

例 用Gauss-Seidel ||x(k+1)-x (k)||2

-1 4 -1 0 -1 0 x2 5

0 -1 4 -1 0 -1 x3 = -2

-1 0 -1 4 -1 0 x4 5

0 -1 0 -1 4 -1 x5 -2

0 0 -1 0 -1 4 x6 6

四. 实验代码：

function[x,k,flag]=Gau_Seid(A,b,delta,max1)

%Gauss-Seidel迭代法

%A为方程组的系数矩阵

%b为方程组的右端

%delat为精度要求，实验中为0.0001

%max1为最大迭代次数，实验中设为50

%x为方程组的解;%k为迭代次数;

%flag为指标变量 flag='ok！' 表示迭代收敛到指标要求 ；flag='fail'表示迭代失败 if nargin

if nargin

n=length(A);k=0;

x=zeros(n,1);y=zeros(n,1);flag='ok';

while 1

y=x;

for i=1;n

z=b(i);

for j=1:n

if j~=i

z=z-A(i,j)*x(j);

end

end

if (abs(A(i,j))

flag='fail';return;

end

z=z/A(i,i);x(i)=z;

end

if norm(y-x,inf)

break;

end

k=k+1;

end

五. 实验结果:

>> A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0 ;0 -1 4 -1 0 -1 ;

-1 0 -1 4 -1 0;0 -1 0 -1 4 -1; 0 0 -1 0 -1 4]

A =

4 -1 0 -1 0 0

-1 4 -1 0 -1 0

0 -1 4 -1 0 -1

-1 0 -1 4 -1 0

0 -1 0 -1 4 -1

0 0 -1 0 -1 4

>> b=[0 5 -2 5 -2 6]

b =0 5 -2 5 -2 6

>> delta=0.0001

delta =1.0000e-004

K>> Gau_Seid(A,b,delta,max1)

11 if nargin

K>> Gau_Seid(A,b,delta,max1)

11 if nargin

六. 实验总结:

通过这次实验, 可分析如下结果:

1) 优点分析:

1. 程序容错性好，对于能够求解和不能够求解的方程程序都会进行处理并给予回复.

2. 程序简单易懂，可读性强。

3. 执行效率高，即使是大型矩阵也可以很快求解

2)缺点分析:

1. 程序每次只能对一个方程组进行求解，对于矩阵A ，b ，程序能一次调用就全部求解。

2. 算法没能很好的利用MATLAB 强大的矩阵运算功能来减少代码，使得代码冗长。

实验四 插值法

————三次样条插值法

一. 实验目的:

1． 学会用三次样条函数作插值计算

2． 通过应用和实际编程计算进一步了解该方法的原理和优缺点。

二. 实验算法:

实验算法（求三次样条插值勤函数S （x ）在一列点x 1,x 2, „x m 上的近似值）；

（1） 输入x i ,y i ,(i=0,1,„,n), 边界条件y 0,y n , 及x j (j=1,2,„m).

（2） 计算h i .h i =xi ,x i-1 i=1,2,„n.

（3） 按下列公式计算μi , λi,g i , i=1,2,„n-1

μi =hi /(hi +hi+1), λi , μi

g i =6/(hi +hi+1)(yi+1-y i )/hi+1-(yi -y i-1)/hi )

(4) 按下列公式计算g 0,g n ,

g 0=(6/h1)((y1-y 0)/h1-y 0`)

g n =6/hn (yn `-(yn -y n-1/hn ))

(5) 用追赶法求解如下方程组得M i , i=0,1,„n,

0 0

μ1 2 λ1 11

„ „ „﹕

μn-1 2 λg n-1

n n

（6）确定点x j (j=1,2,„m) 所在的区间[xi-1,xi],并按下列公式计算S （xj ）; S （x j ）3322=Mi-1(xi -x) +Mi ((x-xi-1) /6hi )+(yi-1-(Mi-1/6)hi )((xi -x)/hi )+(yi -(Mi /6)hi )((x-xi-1)/hi ) (x∈[xi-1,xi],i=1,2,„n)

三. 实验内容:

已知直升飞机旋转机翼外开曲线轮廓线（如图）上的某些弄值点（见表）及端点处的一阶导数数值

y`(x0)=1.86548, y(x18)=-0.046115

试计算该曲线上横坐标为

2，4，6，12，16，30，60，110，180，280，280，400，515

四. 实验代码:

function S=csfit(x,y,dx0,dxn)

%x,y分别为n 个接点的横坐标所组成的向量及纵坐标组成的向量

%dx(),dn分别为S 的导数在x0,xn 处的值

n=length(x)-1;

h=diff(x);

d=diff(y)./h;

a=h(2:n-1);

b=2*(h(1:n-1)+h(2:n));

c=h(2:n);

u=6*diff(d);

b(1)=b(1)-h(1)/2;

u(1)=u(1)-3*(d(1)-dx0);

for k=2:n-1

temp=a(k-1)/b(k-1);

b(k)=b(k)-temp*c(k-1);

u(k)=u(k)-temp*u(k-1);

end

m(n)=u(n-1)/b(n-1);

for k=n-2:-1:1

m(k+1)=(u(k)-c(k)*m(k+2))/b(k);

end

m(1)=3*(d(1)-dx0)/h(1)-m(2)/2;

m(n+1)=3*(dxn-d(n))/h(n)-m(n)/2;

for k=0:n-1

S(k+1,1)=(m(k+2)-m(k+1))/(6*h(k+1));

S(k+1,2)=m(k+1)/2;

S(k+1,3)=d(k+1)-h(k+1)*(2*m(k+1)+m(k+2))/6;

S(k+1,4)=y(k+1);

End

五. 实验结果

>> x=[0,1,2,3];y=[0,0.5,2,1.5];dx0=0.2;dxn=-1;

>> s=csfit(x,y,dx0,dxn)

s =

0.4731 -0.1731 0.2000 0

-1.0192 1.2462 1.2731 0.5000

0.6558 -1.8115 0.6558 2.0000

六. 实验总结:

通过本次实验, 考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，L n (x ) 是否也更加靠近被逼近的函数。龙格（Runge ）给出一个例子是极著名并富有启发性的。项式插值是不收敛的，即插值的节点多，效果不一定就好。

按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数, 可以用MATLAB 的函数“spline ”作此函数的三次样条插值。

实验五 数值积分法

一. 实验目的:

1. 学会用复合辛普森公式计算某积分的近似值；

2. 学会用半步长复合梯形法计算某积分的近似值；

3. 通过实际编程计算进一步了解方法的原理和优缺点。

二. 实验算法:

实验算法：

1. 复合辛普森公式

∫a f（x ）≈h/6[f(a)+4∑f(Xk+1/2)+2∑f(Xk )+f(b)]

其中X k =a+k*h,(k=0，1，...,n),h=（b-a ）/n,Xk+1/2=Xk +h/2。

2. 半步长复合梯形法的递推算法：

T1=（b-a ）/2*[f(a)+f(b)]

T2=1/2*Tn+h/2∑f(Xk+1/2)

其中X k =a+k*h, (k=0，1，...,n-1), h=（b-a ）/n,Xk+1/2=Xk +h/2。 b

三. 实验内容：

1. 用复合辛普森公式计算积分∫0 48（1+cos2x ）1/2 dx ，是误差不超过10-4

2. 用半步长复合梯形法的递推算法计算积分

∫01（1-e -x ）/x dx

使误差不超过10-4

四. 实验代码:

function s=trapr1(f,a,b,n)

h=(b-a)/n;

s=0;

for k=1:n-1

x=a+h*k;

s=s+feval('lp',x);

end

s=h*(feval('lp',a)+feval('lp',b))/2+h*s;

function z=lp(x)

z=1/(1+x^2);

function z=lp(x)

z=(1-exp(-x)^(1/2))/x;

五. 实验结果:

>> trapr1('lp',0.00000001,1,10)

ans =

0.4439

>>

>> trapr1('lp',-1,1,10)

ans =

1.5675

六. 实验总结:

通过本次实验可知, 线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。

高维空间中的积分，如果维数不很高且积分区域是规则的或者能等价地写成多重积分的形式，可以用一元函数积分的数值方法来计算高维空间的积分。蒙特卡罗方法对计算复杂区域甚至不连通的区域上的积分并没有特殊的困难。

实验六 常微分方程的数值解法

一. 实验目的:

1. 学会用改进的欧拉法求一阶常微分方程初值问题的数值解；

2. 学会用经典R-K 法求一阶常微分方程初值问题的数值解；

3. 通过编程计算进一步了解方法的优缺点

二. 实验算法:

实验算法：

1. 改进的欧拉公式：

T 1=y i +h*f(x i ,y i )

T 2=y i +h*f(x i +1,T1)

y i+1=(T1+T2)/2

2. 经典R-K 公式：

y i+1=y i +h/6*(K1+2*K2+2*K3+K4)

K 1=f(x i ,y i )

K 2=f(x i +h/2,y i +h/2*K1)

K 3=f(x i +h/2,y i +h/2*K2)

K 4=f(x i +h,y i+h*K3)

三. 实验内容：

用下列方法求

-2 dy/dx=2/3*x*y,x ∈[0,1]

y(0)=1

21/3的数值解（取h=0.1），并将计算结果与准确解y=(1+x) 进行比较：

（1） 改进的欧拉法

（2） 经典R-K 法

四. 实验代码:

function E=Euler(fun,x0,y0,xN,N)

x=zeros(1,N+1);y=zeros(1,N+1);

x(1)=x0;y(1)=y0;

h=(xN-x0)/N;

for n=1:N

x(n+1)=x(n)+h;

z0=y(n)+h*feval(fun,x(n),y(n));

y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0));

end

T=[x',y']

function z=f(x,y)

z=x+y;

function R=rk4(f,a,b,ya,N)

h=(b-a)/N;

T=zeros(1,N+1);

Y=zeros(1,N+1);

T=a:h:b;

Y(1)=ya;

for j=1:N

k1=h*feval(f,T(j),Y(j));

k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);

k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);

k4=h*feval(f,T(j)+h,Y(j)+k3);

Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;

end

R=[T' Y']

function rk=rk(x,y)

rk=y-2*x/y;

五. 实验结果:

%取迭代次数为5

>> Euler('f',0,1,0.1,5)

T =

0 1.0000

0.0200 1.0204

0.0400 1.0416

0.0600 1.0637

0.0800 1.0866

0.1000 1.1103

function a=b(x,y)

a=(2/3)*x/y^2;

>> Euler('b',0,1,0.1,5)

T =

0 1.0000

0.0200 1.0001

0.0400 1.0005

0.0600 1.0012

0.0800 1.0021

0.1000 1.0033

%R-K法

%迭代次数为10

>> rk4('rk',0,1,1,10)

R =

0 1.0000

0.1000 1.0954

0.2000 1.1832

0.3000 1.2649

0.4000 1.3416

0.5000 1.4142

0.6000 1.4832

0.7000 1.5492

0.8000 1.6125

0.9000 1.6733

1.0000 1.7321

ans =

0 1.0000

0.1000 1.0954

0.2000 1.1832

0.3000 1.2649

0.4000 1.3416

0.5000 1.4142

0.6000 1.4832

0.7000 1.5492

0.8000 1.6125

0.9000 1.6733

1.0000 1.7321

实验

>> rk4('b',0,1,1,10)

R =

0 1.0000

0.1000 1.0033

0.2000 1.0132

0.3000 1.0291

0.4000 1.0507

0.5000 1.0772

0.6000 1.1079

0.7000 1.1422

0.8000 1.1793

0.9000 1.2187

1.0000 1.2599

ans =

0 1.0000

0.1000 1.0033

0.2000 1.0132

0.3000 1.0291

0.4000 1.0507

0.5000 1.0772

0.6000 1.1079

0.7000 1.1422

0.8000 1.1793

0.9000 1.2187

1.0000 1.2599

六. 实验总结:

通过本次实验, 可得到如下结果：

[结果分析]：

从下面的图来看，我们虽然用的是不同的方法，但是求出来的结果却是非常吻合。

[程序的优缺点]：

优点：函数通用性好（特别是经典4阶的Runge-Kutta 方法），只需要稍作修改就可以直接用到其他的微分方程组上。

缺点：简单Euler 方程精度不是很高。

改进的欧拉法也称为梯形法。不难发现，欧拉公式是关于y n +1的显式，即只要已知y n , 经过一次计算便可得y n +1的值，而改进的 欧拉公式是以y n +1的隐式方程给出，不能直接得到y n +1。隐式方程通常用迭代法求解，而迭代过程的实质是逐步显式化。

R-K 方法不是通过求导数的方法构造近似公式，而是通过计算 不同点上的函数值，并对这些函数值作线性组合，构造近似公式，再把近似公式与解的泰勒展开 式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。