数学与信息工程学院
实 验 报 告
课程名称: 数值分析
实 验 室:
实验台号:
班 级:
姓 名:
实验日期: 2012 年 3 月 14 日
第二篇:数值分析实验报告册
实验名称: Lagrange插值(实验一)
实验目的:
掌握Lagrange插值数值算法,能够根据给定的函数值表达求出插值多项式和函数在某一点的近似值。
实验准备:
1. 在开始本实验之前,请回顾教科书的相关内容;
2. 需要一台准备安装Windows XP Professional操作系统和装有数学软件的计算机。
实验内容及要求
已知数据如下:
要求:
试用Lagrange插值多项式求时的函数近似值.
实验过程:
编写Matlab函数M文件Lagrange如下:
function yy=lagrange(x,y,xi)
m=length(x); n=length(y);
if m~=n,error('向量x与y的长度必须一致');end
for k=1:length(xi)
s=0;
for i=1:m
z=1;
for j=1:n
if j~=i
z=z*(xi(k)-x(j))/(x(i)-x(j));
end
end
s=s+z*y(i);
end
yy=s
end
在命令窗口调用函数M文件lagrange,输出结果如下:
>>x=[0.56160, 0.56280, 0.56401, 0.56521];
>>y=[0.82741, 0.82659, 0.82577, 0.82495];
>>xi=[0.5626, 0.5635, 0.5645];
>>yi= lagrange (x,y,xi)
yi=
0.8628 0.8261 0.8254
实验总结(由学生填写):
教师对本次实验的评价(下面的表格由教师填写):
实验名称: 曲线拟合的最小二乘方法(实验二)
实验目的:
掌握最小二乘方法,并能根据给定数据求其最小二乘一次或二次多项式,然后进行曲线拟合。
实验准备:
1. 在开始本实验之前,请回顾教科书的相关内容;
2. 需要一台准备安装Windows XP Professional操作系统和装有VC++6.0的计算机。
实验内容及要求
炼钢是个氧化脱碳的过程,钢液含碳量的多少直接影响冶炼时间的长短,下表是某平炉的生产记录,表中为实验次数,为全部炉料熔化完毕时钢液的含碳量,为熔化完毕至出钢所需的冶炼时间(以分为单位).
将所数据通过图示方法绘在坐标纸上,观察数据点的分布情况,然后进行曲线拟合.
实验过程:
编写Matlab命令文件如下:
function t=zxecnh(x,y)
x=[-1.00 -0.75 -0.50 -0.25 0];
y=[-0.2209 0.3295 0.8826 1.4392 2.0003];
syms sumx sumy sumxy sumx2;
sumx=0
sumy=0
sumxy=0
sumx2=0
for i=1:5
sumx=sumx+x(1,i);
sumy=sumy+y(1,i);
sumxy=sumxy+(x(1,i))*(y(1,i));
sumx2=sumx2+(x(1,i))*(x(1,i));
end
sumx
sumy
sumxy
sumx2
b=(sumy*sumx-sumxy*5)/(sumx*sumx-sumx2*5)
a=(sumy-b*sumx)/5
scatter(x,y)
hold on
plot(x,a+b*x)
在命令窗口调用函数M文件zxecnh,输出结果如下:
b = 2.2208
a = 1.9966
实验总结(由学生填写):
教师对本次实验的评价(下面的表格由教师填写):
实验名称: Romberg积分法(实验三)
实验目的:
掌握Romberg算法,并能根据给定的精度要求计算定积分。
实验准备:
1. 在开始本实验之前,请回顾教科书的相关内容;
2. 需要一台准备安装Windows XP Professional操作系统和装有VC++6.0的计算机。
实验内容及要求
用Romberg法求函数积分,精度为.
实验过程:
编写函数M文件Romberg如下:
funciton 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)
>>format long;Romberg(inline( 'sin(x)./x'),eps,1,0.5e-6);format short;
T=
t=0.94608307038722
实验总结(由学生填写):
教师对本次实验的评价(下面的表格由教师填写):
实验名称: 常微分方程的差分方法(一、二)(实验四)
实验目的:
掌握欧拉方法与改进的欧拉方法、三阶、四阶龙格—库塔方法,能够用三阶和四阶经典公式求解微分方程。
实验准备:
1. 在开始本实验之前,请回顾教科书的相关内容;
2. 需要一台准备安装Windows XP Professional操作系统和装有VC++6.0的计算机。
实验内容及要求
(1)求初值问题:
.
(2)取步长,写出用经典四阶龙格—库塔方法求解初值问题
的计算公式.
实验过程:
(1)编写函数M文件改进的欧拉方法如下:
function E=euler(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;
T=zeros(1,M+1);
Y=zeros(1,M+1);
Yp=0;
Yc=0;
T=a:h:b;
Y(1)=ya;
for j=1:M
Yp=Y(j)+h*(T(j)-Y(j))/2;
Yc=Y(j)+h*(T(j)-Yp)/2;
Y(j+1)=(Yp+Yc)/2;
end
T
Y
在命令窗口调用函数M文件euler(0,1,1,10),输出结果:
T =0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
Y = 1.0000 0.9512 0.9098 0.8752 0.8471 0.8253 0.8095 0.7992 0.7944 0.7947 0.7998
(2)编写函数M文件4阶龙格-库塔方法如下:
function R=rk4(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;
T=zeros(1,M+1);
Y=zeros(1,M+1);
T=a:h:b;
Y(1)=ya;
for j=1:M
k1=h*(T(j)*(cos(T(j)+Y(j))));
k2=h*((T(j)+h/2)*(cos((T(j)+h/2)+(Y(j)+k1/2))));
k3=h*((T(j)+h/2)*(cos((T(j)+h/2)+(Y(j)+k2/2))));
k4=h*((T(j)+h)*(cos((T(j)+h)+(Y(j)+k3))));
Y(j+1)=Y(j)+(k1+k2+k3+k4)/6;
end
T
Y
在命令窗口调用函数M文件rk4 (1,4,0,10),输出结果:
T =1.0000 1.3000 1.6000 .9000 2.2000 2.5000 2.8000 3.1000 3.4000 3.7000 4.0000
Y =0 0.0762 0.0825 0.0071 -0.1447 -0.3596 -0.6225 -0.9187 -1.2357 -1.5635 -1.8949
实验总结(由学生填写):
教师对本次实验的评价(下面的表格由教师填写):
实验名称: Newton方法(实验五)
实验目的:
掌握Newton迭代算法,能够根据所给方程求出在某一点附近的根。
实验准备:
1. 在开始本实验之前,请回顾教科书的相关内容;
2. 需要一台准备安装Windows XP Professional操作系统和装有VC++6.0的计算机。
实验内容及要求
公元1225年,Lenardo宣布他求得方程
的一个根.当时颇为轰动,但无人知道他是用什么方法得到的.现在,请你试试用Newton迭代法求解这个结果.
实验过程:
编写Matlab函数M文件Newton如下:
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=le-4;end
X=x0;x0=x+2﹡e;k=0;
Fprintf(‘It.no=%2d x%[2d]=%12.9f\n`,k,k,x)
While abs(x0-x)>e£k<N
K=K+1
X0=x;x=x0-feval(fname,x0)/feval(dfname,x0);
Fprintf(It.no=%2d x[%2d]=%12.9f\n`,k,k,x)
End
If k==N,fprintf(`已达到迭代次数上限`);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);
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
五、弦截法
例4.15 用弦截法求方程f(x)=x(x+1)2-1=0在0.4附近的一个实根.初始值x0=0.4,x1=0.6,精确指4位有效数字.
编写Matlab函数M文件XianjieMethod和chap4-fun如下:
Function f=XianjieMethod(x0,x1)
X2=x1-chap4-fun(x1)*(x1-x0)/chap4—fun(x1)—chap—fun(x0));
Eps=le—4;
N==0;
rprintf('迭代次数 x—n feval(x—n)\n'')
rprintf('n=%3.0f x--%d=%10.5f %10.6e\n',n,n,xx0,chap4—fun(x0))
while abs(x1—x0)>eps&(n<600)
x0=x1;x1=x2;
x2=x1—chap4—fun(x1)*(x1-x0)/(chap4—fun(x1)—chap4-fun(x0));
n=n+1
fprintf('n=%3.0f x--%10.5f %10.6e\n',n,n,x0,chap4—fun(x0))
end
fprintf('\n迭代次数n=%3.0f x﹡=%10.5f',nx0)
function f=chap4—fun(x)
f=x*(x+1)^2-1;
取x0=0.4, x1=0.6在命令窗口调用函数M文件XianjieMeethod输出结果如下:
>>XianjieMethod(0.4,0.6)
迭代次数 x-n feval(x-n)
n=0 x-0=0.40000 -2.160000e-001
n=1 x-1=0.60000 5.360000e-001
n=2 x-2=0.45745 -2.831381e-002
n=3 x-3=0.46460 -3.410914e-003
n=4 x-4=0.465558 2.698819e-005
迭代次数n=4 x*=0.46558
实验总结(由学生填写):
教师对本次实验的评价(下面的表格由教师填写):
实验名称: 雅可比迭代方法(实验六)
实验目的:
掌握列主元的高斯消去法思想,能够利用列主元的高斯消去法求解任意阶数的线性方程组。
实验准备:
1. 在开始本实验之前,请回顾教科书的相关内容;
2. 需要一台准备安装Windows XP Professional操作系统和装有VC++6.0的计算机。
实验内容及要求
用的雅可比迭代方法求解下列方程组:
.
实验过程:
编写Matlab函数M文件jacobi如下:
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<delta)|(relerr<delta)
break
end
end
X=X';
在命令窗口调用函数M文件jacobi([10 0 1 -5;1 8 -3 0;3 2 -8 1;1 -2 2 7],[-7 11 23 17]',[0 0 0 0]',1.0e-005,10000)输出结果如下:
ans = 1.0000 0.5000 -2.0000 3.0000
实验总结(由学生填写):
教师对本次实验的评价(下面的表格由教师填写):
实验名称: 列主元的高斯消去法(实验七)
实验目的:
掌握列主元的高斯消去法思想,能够利用列主元的高斯消去法求解任意阶数的线性方程组。
实验准备:
1. 在开始本实验之前,请回顾教科书的相关内容;
2. 需要一台准备安装Windows XP Professional操作系统和装有VC++6.0的计算机。
实验内容及要求
解下列方程组
实验过程:
编写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
实验总结(由学生填写):
教师对本次实验的评价(下面的表格由教师填写):