数值分析实验报告

数值分析 实验报告

系 别 专 业 班级学号 姓 名 指导教师

电子信息系 计算机科学与技术

2011年 6 月20

日

数值分析实验报告

实验一

实验题目：

编写一个拉格朗日插值函数，对不多于9个点的插值节点都可以求出插值函数，任意给定输入x值都可以求出y值。 例如：（0,1）,(1,1) ,(4,2) ,(9,3) ,(16,4) ,(25,5) ,(36,6) ,(49,7) ,(64,8)，这实际上是一个平方根函数随意输入x在[0，64]上的值都可以输出一个y值。 进一步可以在这区间上画出图像。

程序代码：

function yi = Lagrange(x,y,xx) %Lagrange 差值多项式，其中 %x为向量，全部差值节点， %y为向量，差值节点处的函数值， %xi为标量，为被估计函数的自变量， %yi为xi处的函数估计值，

n = length(x);m = length(y); %输入差值节点与它的函数值应有相同个数， if n ~=m

error('The lengths of X and Y must be equal!'); return; end

p = zeros(1, n); for k =1 : n t = ones(1, n); for j = 1 : n if j ~=k

%输入差值节点必须互异

if abs(x(k)-x(j)) < eps error('the DATA is error!'); return; end

t(j) = (xx - x(j)) / (x(k) - x(j)); end end

p(k) = prod(t); end

yi = sum(y .* p); end

测试数据：

数值分析实验报告

>> clear

>> x=[0 ; 1 ; 4 ; 9 ; 16 ; 25 ; 36 ; 49 ; 64]; >> y=0:8;

>> f=interp1(x,y,y); >> F=interp1(x,y,y,'spline'); >> f,F f =

0 1.0000 1.3333 1.6667 2.0000 2.2000 2.4000 2.6000 2.8000 F =

0 1.0000 1.5601 1.8402 2.0000 2.1706 2.3682 2.5806 2.7953 >> clear

>> x=[0 ; 1 ; 4 ; 9 ; 16 ; 25 ; 36 ; 49 ; 64]; >> y=0:8; >> t=0:0.1:70;

>> f=interp1(x,y,t,'spline'); >> plot(x,y,'k+',t,f,x,y,'r') >>xlabel('x') >>ylabel('y') >>titel('y=x^1/2') 实验结果：

数值分析实验报告

f =

0 1.0000 1.3333 1.6667 2.0000 2.2000 2.4000 2.6000 2.8000 F =

0 1.0000 1.5601 1.8402 2.0000 2.1706 2.3682 2.5806 2.7953

实验二

实验题目：用不同的方法计算积分 0xlnxdx

1

4 9

取不同的步长h，分别用复合梯形公式及复合辛普森公式计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善。

程序代码：

数值分析实验报告

1.编写一个函数实现被积函数 Shiyan2.m

function [ y ] = shiyan2( x )

%UNTITLED Summary of this function goes here % Detailed explanation goes here y=log(x).*sqrt(x); end

2．使用matlab画出函数曲线 在matlab中输入

>> x = 0.01:0.001:1 >> y = shiyan2(x) >> plot(x,y,'-b')

3．编写复合梯形公式计算积分的函数 function [ y ] = echelon( h )

%UNTITLED3 Summary of this function goes here % Detailed explanation goes here

数值分析实验报告

n = length(h); y=[0,0];

for i = 1:n; temp=0.001; result = 0;

while temp+h(i)<1

x=(shiyan2(temp)+shiyan2(temp+h(i)))*h(i)/2; result=result+x; temp = temp+h(i); end

result = result +(shiyan2(temp)+shiyan2(1))*h(i)/2; y(i) = result; end end

输入一组步长

h

；

调用符合梯形积分函数计算积分结果

计算积分结果的误差

数值分析实验报告

5.编写使用辛普森公式计算积分的函数 function [ y ] = fsimp(h)

%UNTITLED5 Summary of this function goes here % Detailed explanation goes here temp=0.001; result=0;

while temp+h<1

x=4*shiyan2(temp+h/2)+2*shiyan2(temp+h); result=result+x; temp=temp+h; end

result=[shiyan2(temp)+shiyan2(1)+result]*h/6; y=result; end

trapz(x,y),使用梯形公式计算积分，参数x为定义域内一系列值构成的向量，y=f(x) 例如：

>> x=0.1:0.1:1; >> y=shiyan2(x); >> trapz(x,y)

数值分析实验报告

ans =

-0.3807 >>

quad(fun,a,b)使用辛普森公式计算积分 fun代表函数句柄，a,b为积分上下限 >> s=@shiyan2 s =

@shiyan2

>> quad(s,0,1)

ans =

-0.4444

>>plot(y,’DisplayName’,’y’,’TDatasource’,’y’

);figure(gcf)

实验三

数值分析实验报告

实验题目：

用LU分解和列主元消去法解线性方程组

7 10

32.099999 5 1

1 2

1 x1 8

x 5.900001 62 2

5 1 x3 5

02 x4 1 0

输出Ax=b中系数A=LU分解的矩阵L和U，解向量x及detA；列主元法的行交换

次序，解向量x及detA；比较两种方法所得的结果。

程序代码：

1.新建函数lux.m

function [L,U,x]=lux(A,b) [n,n]=size(A);

p=eye(n);//eye函数，生成一个单位矩阵 for k=1:n-1

[r,m]=max(abs(A(k:n,k))); m=m+k-1; if(A(m,k)~=0) if(m~=k)

A([k m],:)=A([m k],:); p([k m])=p([m k]); end

for i=k+1:n

A(i,k)=A(i,k)/A(k,k); j=k+1:n;

A(i,j)=A(i,j)-A(i,k)*A(k,j); end end end

L=tril(A,-1)+eye(n,n);提取下三角矩阵 U=triu(A);提取上三角矩阵 newb=p*b; y=zeros(n,1);全0 for k=1:n j=1:k-1;

y(k)=(newb(k)-L(k,j)*y(j))/L(k,k); end

x=zeros(n,1); for k=n:-1:1 j=k+1:n;

x(k)=(y(k)-U(k,j)*x(j))/U(k,k); end

数值分析实验报告

输入矩阵A

>> A=[10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2] 输入b

B=[8;5.900001;5;1] 调用函数lux计算结果 >> [l,u,x]=lux(A,B) 调用函数det求行列式

>> det(A)

数值分析实验报告

二 高斯列主元消去法

function x=gauss_lie(A,b) %采用高斯列主元法求解方程组Ax=b n=length(b); p=1:n;lu=A; y=[]; for k=1:n

[c,i]=max(abs(lu(k:n,k))); ik=i+k-1; if ik~=k

m=p(k);p(k)=p(ik);p(ik)=m;

ck=lu(k,:);lu(k,:)=lu(ik,:);lu(ik,:)=ck; end if k==n break; end

lu(k+1:n,k)=lu(k+1:n,k)/lu(k,k);

lu(k+1:n,k+1:n)=lu(k+1:n,k+1:n)-lu(k+1:n,k)*lu(k,k+1:n); end

l=diag(ones(n,1))+tril(lu,-1); u=triu(lu);

y(1)=b(p(1)); for i=2:n

y(i)=b(p(i))-l(i,1:i-1)*y(1:i-1)'; end

x(n)=y(n)/u(n,n); for i=n-1:-1:1

x(i)=(y(i)-u(i,i+1:n)*x(i+1:n)')/u(i,i); end x=x'; end

测试数据：A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]

数值分析实验报告

b=[8 5.900001 5 1]

实验结果：

> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2] A =

10.0000 -7.0000 0 1.0000 -3.0000 2.1000 6.0000 2.0000 5.0000 -1.0000 5.0000 -1.0000 2.0000 1.0000 0 2.0000

>> b=[8 5.900001 5 1] b =

8.0000 5.9000 5.0000 1.0000 >>

> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];[L,U,index]=LU_Decoml(A) L =

1.0e+006 *

0.0000 0 0 0 -0.0000 0.0000 0 0 0.0000 -2.5000 0.0000 0 0.0000 -2.4000 0.0000 0.0000 U =

1.0e+007 *

0.0000 -0.0000 0 0.0000 0 -0.0000 0.0000 0.0000 0 0 1.5000 0.5750 0 0 0 0.0000

数值分析实验报告

index =

1

实验四

实验题目：

下列方程的实根： （1）x（2）x

2

3x 2 ex 0 2x2 10x 20 0

3

要求：（1）设计一种不动点迭代法，要使迭代序列收敛，然后再用Atiken加速迭代，计算到

xk xk 1 10 8

为止；（2）用牛顿迭代法，同样计算到

xk xk 1 10 8

，输出迭代初值和迭代次

数k，比较方法的优劣。

数值分析实验报告

程序代码：

1． 构造迭代公式：y=(x2+2-ex)/3

function [ y ] = shiyan4_1( x )

%UNTITLED3 Summary of this function goes here % Detailed explanation goes here y=(x*x+2-exp(x))/3; end

2． 参考以下不动迭代法计算非线性方程的代码 function [ result ] = Iteration( x )

%UNTITLED Summary of this function goes here % Detailed explanation goes here x1=shiyan4_1(x); x0 = x;

while(abs(x1-x0)>1e-8) temp = shiyan4_1(x1); x0=x1; x1=temp; end

result = x1; end

3． 参考以下使用艾特金加速迭代计算线性方程的代码

function [ result ] = steffensen( x )

%UNTITLED2 Summary of this function goes here % Detailed explanation goes here x1=shiyan4_1(x); x0 = x;

while(abs(x1-x0)>1e-8) temp = sub(x1); x0=x1; x1=temp; end

result = x1;

function b=sub(a)

b=a-sqrt(shiyan4_1(x)-x)./(shiyan4_1(shiyan4_1(x))-2*shiyan4_1(x)+x); end

数值分析实验报告

end

测试数据：

>>X=0.5;

>>Y=Iteration( x ); 结果： y= 0.2575 测试数据：

>> l=Newton(0.5) 结果：l =

0.2575