⒈目的

数值计算的应用范围十分广泛，作为用计算机解决实际问题的基础、桥梁和纽带，是架设在数学与计算机之间的一条通道．数值计算方法实验将帮助学生理解并掌握数值计算方法这门专业课的基本理论和基本内容，准确地将所求解的数学模型简化成一系列算术或逻辑运算，上机实现其数值求解．通过实验，加深学生对一些重要算法的理解，提高学生的编程能力与解决实际问题的能力，培养学生应用计算方法解决工程计算的能力，以期达到初步的科学计算和研究的目的．

⒉实验任务分解

通过一些实例掌握数值计算的基本计算方法，针对基本计算方法，运用数值计算软件Matlab进行数学编程求解，并对问题中的各有关变量进行分析、调试和推广．

⒊实验环境介绍

计算机房

⒋实验时数

16学时

实验一

⒈实验目的与要求

通过对具体实例，理解并掌握Lagrange插值的数值算法，能够根据实例中给定的函数值表求出插值多项式和函数在某一点的近似值，并运用Matlab编程求解．

⒉实验内容

已知数据如下：

	10	11	12	13
	2.3026	2.3979	2.4849	2.5649

要求：

试用分别用线性插值、二次插值、三次插值分别计算在点 时的函数近似值．

⒊思考题

（1）       Lagrange插值法在编程中最关键的问题是什么？

（2）       编写一个用Newton前插值公式计算函数值的程序，要求先输出差分表，再计算 点的函数值，并应用于下面的问题：

	20	21	22	23	24
	1.30103	1.32222	1.34242	1.36173	1.38021

求 时的三次插值多项式的值．

（3）       在某海域测得一些点 处的水深 见如下数据，在矩形

区域内画出海地面的图形．

 

实验二

⒈实验目的与要求

掌握最小二乘方法，并能根据给定数据求其最小二乘一次或二次多项式，然后进行曲线拟合。

⒉实验内容

已知一组实验数据如下：．

	1	2	3	4
	2	4	6	8
	2	11	28	40

试用最小二乘法求一次拟合多项式，并将拟合曲线画出来．

⒊思考题

混凝土的抗压强度随养护时间的延长而增加，现将一批混凝土作成12个试块，记录了养护日期x（日）及抗压强度y（kg/cm2）的数据： 

已知建立的非线性回归模型为y=a+k1*exp(m*x)+k2*exp(-m*x)，要求对得到的模型和系数进行检验。 

注明：此题中的+r代表加上一个[-0.5,0.5]之间的随机数 

 

实验三

⒈实验目的与要求

通过对具体实例中给出了积分问题，学会运用Romberg算法在给定精度下计算定积分问题，并运用Matlab编程求解．

⒉实验内容

用Romberg法求函数积分 ，精度为 ．

⒊思考题

（1）用Matlab求解广义积分：

     （2）用复合梯形公式计算下面积分，取 ：

并给出数值积分结果与精度 之间的误差．

实验四

⒈实验目的与要求

掌握改进的欧拉方法、三阶、四阶方法，会用三阶和四阶经典的龙格—库塔格式求解所给实例的初值问题，并运用Matlab编程求解．

⒉实验内容

（1）用改进的欧拉方法求初值问题：

       ．

（2）取步长 ，写出用经典四阶龙格—库塔方法求解初值问题

的计算公式．

⒊思考题          

（1）用欧拉算法编程求解解一阶微分方程初值问题：

（2）小船渡过宽为 的河流A-B， 为河水流速， 为静水速度，渡船时船头始终指向B点；

a）建立小船的航线方程求其解析解；

b）设 m， ， ，用数值解法求渡河所需要时间和任意时刻小船及航行的曲线，作图，并与解析解比较．若流速

，其结果如何？

实验五

⒈实验目的与要求

通过求解方程，学会运用Newton迭代法求解方程，并运用Matlab编程求解．

⒉实验内容

公元1225年，Lenardo宣布他求得方程

的一个根 ．当时颇为轰动，但无人知道他是用什么方法得到的．现在，请你试试用Newton迭代法求解这个结果．

⒊思考题          

(1)        用弦截法求 在0.4附近的一个实根，初始值 ，精确到4位有效数字．

(2)        用对分区间法求解方程

在区间 上的误差小于 的一个根，并记录对分区间的次数．

实验六

⒈实验目的与要求

掌雅可比迭代方法和高斯—赛德尔迭代方法思想，学会用该方法求解方程组，并比较收敛速度。

⒉实验内容

用雅可比迭代方法和高斯—赛德尔迭代方法求解下列方程组：

，

并比较收敛速度。

⒊思考题          

    用超松弛法求解上述方程组，并与高斯—赛德尔迭代方法的速度进行比较。

实验七

⒈实验目的与要求

掌握列主元的高斯消去法思想，能够利用列主元的高斯消去法求解任意阶数的线性方程组。

⒉实验内容

用列主元的高斯消去法解下列方程组

.

 

 

 

 

 

 

 

供参考的实验过程：

实验一

function  y=lagrange(x0,y0,x)

w=length(x0);

n=w-1; 

m=length(x);

for i=1:m

    z=x(i);

    s=0;

        for k=1:n+1

            p=1;

                for j=1:n+1

                    if j~=k

                            p=p*(x(i)-x0(j))/ (x0(k)-x0(j));

                    end

                end

                s=p*y0(k)+s;

        end

        y(i)=s;

end

在命令窗口调用函数M文件lagrange([10 11 12 13],[2.3026 2.3979 2.4849 2.5649],11.75)，即可得到线性插值结果。其它插值类似可得。

 

实验二

function z=zuixiaoercheng(x,y)

syms sumx sumy sumxy sumx2;

sumx=0;

sumy=0;

sumxy=0;

sumx2=0;

n=length(x)

 for i=1:n

    sumx=sumx+x(i);

    sumy=sumy+y(i);

    sumxy=sumxy+(x(i))*(y(i));

    sumx2=sumx2+(x(i))*(x(i));

 end

 sumx

 sumy

 sumxy

 sumx2

 b=(sumy*sumx-sumxy*n)/(sumx*sumx-sumx2*n)

 a=(sumy-b*sumx)/5

  plot(x,y,'*')

hold on

plot(x,a+b*x)

在命令窗口调用函数M文件zuixiaoercheng([2 4 6 8],[2 11 28 40])，输出结果如下：

b = 2.2208

a = 1.9966

图略。

实验三

function t=Romberg(fname,a,b,e)%Rom

%Romberg法求函数的积分

%fname是被积函数，a是上限，b是下限，e为精度（默认1e－4）

if nargin<4 ,e="1e-4;

end

i=1; j=1; h=b-a;

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

T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+0.001*h))*h/2;

T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/ (4^j-1);

while abs(T(i+1,i+1)-T(i,i))>e

i=i+1; h=h/2;

   T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+0.001*h))*h/2;

   for j=1:i

       T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);

   end

end

T

t=T(i+1,j+1)

h

format long;Romberg(inline( 'sin(x)./x＇),eps,1,0.5e－6);format short;

T=

0.92073549240395	0	0	0
0.93979328480618	0.94614588227395	0	0
0.94451352166539	0.94608300406367	0.94608693395179	0
0.94569086358270	0.94608331088847	0.94608306935092	0.94608307038722

 

 

 

 

 

t=0.94608307038722

 

实验四

（1）

function E=gaijineuler(a,b,ya,M)

% Input - a and b are the left and right end points

%       -ya is the initial condition y(a)

%       -M is the number of steps

% Output T is the vector of abscissas and Y is the vector of

%        ordinates

h=(b-a)/M;

Y=zeros(1,M+1);

Yp=0;

Yc=0;

T=a:h:b;

Y(1)=ya;

for j=1:M

    Yp=Y(j)+h*(-Y(j))/2;

    Yc=Y(j)+h*(-Yp)/2;

    Y(j+1)=(Yp+Yc)/2;

end

T

Y

在命令窗口调用函数M文件gaijineuler(0,1.5,1,10)，输出结果：
T = 0.0000    0.1500    0.3000    0.4500    0.6000    0.7500    0.9000    1.0500    1.2000 1.3500    1.5000

Y = 1.0000    1.0778    1.1617    1.2521    1.3495    1.4545    1.5677    1.6897    1.8212  1.9629    2.1156

 （2）

编写Matlab函数M文件Lk如下：

function z=Lk(x0,y0,h,N)

X=zeros(2,N);

function  q=f(x,y)

q=x*cos(x+y);

end

for n=1:N

    x1=x0+h;

   k1=f(x0,y0);

    k2=f(x0+h/2,y0+h/2*k1);

    k3=f(x0+h/2,y0+h/2*k2);

    k4=f(x0+h,y0+h*k3);

    y1=y0+h/6*(k1+2*k2+2*k3+k4);

    X(1,n)=x1;

    X(2,n)=y1;

    x0=x1;

    y0=y1;

end

X

end

在命令窗口调用函数M文件Lk(1,0,0.4,7),输出结果如下:

1.40000000000000   1.80000000000000   2.20000000000000

0.12849071567271   0.03987952838616  -0.23311764949180

2.60000000000000   3.00000000000000   3.40000000000000   3.80000000000000

  -0.61730705974470  -1.05090251367142  -1.49621004434439  -1.93718308658337

实验五

function x=Newton(fname,dfname,x0,e,N)

%用途：Newton迭代法解非线性方程f(x)＝0

%fname和dfname分别表示f(x)及其导数函数的M函数句柄或内嵌函数表达式，

%x0为迭代初值，e为精度（默认值le-4），

%x0为返回数值解，并显示计算过程，设置迭代次数上限N以防发散（默认500次）

if nargin<5 ,

    N=500;

end

if nargin<4 ,

    e=l0^(-4);

end

x=x0;

x0=x+2*e;

k=0;

fprintf('It.no=%2d   x%[2d]=%12.9fn',k,k,x)

while abs(x0-x)>e

   k=k+1

   x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);

   fprintf('It.no=%2d  x[%2d]=%12.9fn',k,k,x)

if k==N,

    fprintf('已达到迭代次数上限')

    break

end

end

在命令窗口编写内嵌函数表达式，并调用函数M文件Newton：

>>fun=inline（’x＾3+2﹡x＾2+10﹡x-20’^｀);

>>dfun=inline( ‘3﹡x＾2+4﹡x+10’);

>>Newton(fun,dfun,1.5,0.5e-6,500);

It.no=0     x[0]=1.500000000

It.no=1     x[1]=1.373626374

It.no=2     x[2]=1.368814820

It.no=3     x[3]=1.368808108

It.no=4     x[4]=1.368808108

 

实验六

function X=jacobi(A,B,P,delta,max1)

%Input -A is an N*N nonsingular matrix

%       -B is an N*1 matrix

%      -P is an N*1 matrix; the initial guess

%      -delta is the tolerance for P

%      -max1 is the maximum number of iterations

%Output -X is an N*1 matrix: the jacobi approximation to the solution of 

%        AX=B

N=length(B);

for k=1:max1

    for j=1:N

        X(j)=(B(j)-(A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N])))/A(j,j);

    end

    err=abs(norm(X'-P));

    relerr=err/(norm(X)+eps);

    P=X';

    if(err

        break

    end

end

X=X';

在命令窗口调用函数M文件jacobi([10 －1 2 0;－1 11 －1 3;2 －1 10 －1;0 3 －1 8],[6 25 －11 15]',[0 0 0 0]',1.0e-005,10000)输出结果如下：

ans =   1.0000    2.0000   -1.0000    1.0000

实验七

编写Matlab命令M文件chap3_1 如下：

Clear

fprintf(‘增广矩阵’)

A=[0.792,0.81,0.9,0.6867;1,1,1,0.8338;1.331,1.21,1,1,1.000]

                                                   %输入增广距阵

%第一次选主元，第三行和第一行交换

fprintf（’第一次选主元后的增广距阵’）

tempo=A(3,:);A(3,:)= A(1,:);A(1,:)=tampo;A

%第一次消元

fprintf(‘第一次消元后的增广距阵’)

A（2,:）= A（2,:）- A（1,:）* A（2,1）/ A（1,1）;

A（3,:）= A（3,:）- A（1,:）* A（3,1）/ A（1,1）;A

%第二次选主元，第三行和第二行交换

fprintf(‘第二次消元后的增广距阵’)

tempo=A(3,:);A(3,:)= A(2,:);A(2,:)=tampo

%第二次消元

fprintf(‘第二次消元后的增广距阵’)

A（3,:）= A（3,:）- A（2,:）* A（3,2）/ A（2,2）;A

%回代求解

fprintf(‘回代求解’)

x(3)=A(3,4)/A(3,3);

x(2)=(A(2.4)-A(2,3)*x(3))/A(2,2);

x(1)=(A(1,4)-A(1,2,3)*x(2:3)’)/A(1,1)

x

在命令窗口调用M文件chap3_1

>> chap3_1

增广距阵

A=

0.7290       0.8100  0.9000  0.6867

1.0000  1.0000  1.0000  0.8338

1.3310  1.2100  1.1000  1.0000

第一次选主元后的增广距阵

A=

1.3310  1.2100  1.1000  1.0000

1.0000  1.0000  1.0000  0.8338

0.7290  0.8100  0.9000  0.6867

第一次消元后的增广距阵

A=

1.3310       1.2100  1.1000  1.0000

0         0.0909  0.1763  0.0825

0         0.1473  0.2975  0.1390

第二次选主元后的增广距阵

A=

1.3310  1.2100  1.1000  1.0000

0  0.1473  0.2975  0.1390

0  0.0909  0.1763  0.0825

第二次消元后的增广距阵

A=

1.3310  1.2100  1.1000  1.0000

 0  0.1473  0.2975  0.1390

           0      0  -0.0101  -0.0033

回带求解

X=

  0.2245  0.2814  0.3279