数值分析报告Word下载.docx
- 文档编号:18697131
- 上传时间:2022-12-31
- 格式:DOCX
- 页数:25
- 大小:57.18KB
数值分析报告Word下载.docx
《数值分析报告Word下载.docx》由会员分享,可在线阅读,更多相关《数值分析报告Word下载.docx(25页珍藏版)》请在冰豆网上搜索。
%err是p0的误差估计
%k是迭代次数
%y=f(p1)
p0,feval('
f'
p0)
fork=1:
N
p1=p0-feval('
p0)/feval('
df'
p0);
err=abs(p1-p0);
p0=p1;
p1,err,k,y=feval('
p1)
if(err<
delta)|(y==0),
break,end
end
方程函数:
functiony=f(x)
y=exp(x)+10*x-2
微分函数:
functiony=df(x)
y=exp(x)+10;
实验结果:
newton('
'
0,5^(-4),10)
p0=0
y=-1
ans=-1
p1=0.0909
err=0.0909
k=1
y=0.0043
p1=0.0905
err=3.8398e-004
k=2
y=8.0727e-008
y=8.0727e-008
ans=0.0905
六.实验总结:
相关MATLAB函数:
x=fzero(fun,x0)返回一元函数fun的一个零点,其中fun为函数句柄,x0为标量时,返回在x0附近的零点;
x0为向量[a,b]时,返回函数在[a,b]中的零点
[x,f,h]=fsolve(fun,x0)返回一元或多元函数x0附近fun的一个零点,其中fun为函数句柄,x0为迭代初值;
f返回fun在x的函数值,应该接近0;
h返回值如果大于0,说明计算结果可靠,否则不可靠
通过本次实验,加深了对非线性问题的迭代法与线性问题迭代法的差别的认识,深入了解了迭代法及初始值与迭代收敛性的关系。
迭代法是求解非线性方程的基本思想方法,与线性方程的情况一样,其构造方法可以有多种多样,但关键是怎样才能使迭代收敛且有较快的收敛速度。
实验二线性代数方程组的解法
--------列主元Gauss消元法
一.实验目的
1.学会用列主元Gauss消去法求解线性代数方程组
2.加深对选取主元素的重要性的认识
3.通过实际变成计算进一步了解该方法的优缺点
二.实验算法:
实验算法(列主元Gauss消去法)
(1)输入数据A和b
(2)对于k=1,2,…..,n-1,做
A,按列选主元。
令|ark|=max|aik|
K<
=i<
=n
如果|ark|=0,则停止计算
B,交换两行。
如果r>
k,则akj<
->
arjj=k,k+1,…..,nbk<
br
C,消元计算。
置mik=aik/akki=k+1,k+2,…,n
aij=aij-mik*akji,j=k+1,k+2,…,n
bi=bi-mik*bki=k+1,k+2,…n
(3)如果ann=0,则停止计算,否则进行回代,对于k=n,n-1,…1,置
xk=(bk-∑akj*xj)/akk(当下标大于上标时,不做求和运算)
(4)输出x1,x2,…,xn,即为方程组Ax=b的解
三.实验内容:
用列主元Gauss消去法求解线性代数方程组
1.1348
3.83261.16513.4017x19.5342
1.53011.78752.53301.5435x26.3941
3.41294.93178.76431.3142x3=18.4231
1.23714.999810.67210.0147x416.9237
4.实验代码及结果:
functionu=gauss(a,b)
n=length(b);
%success=1;
%判断gauss消去法是否成功,1成功,0失败
b_tmp=1:
n;
%记录列变换的情况
%消去过程begin
fork=1:
n-1
%选主元begin
max=abs(a(k,k));
%当前主元的最大值
maxL=k;
%当前主元所在的行
maxC=k;
%当前主元所在的列
forj=k:
n
fori=k:
ifabs(a(j,i))>
max
%[maxL,maxC]=[j,i];
maxL=j;
maxC=i;
max=abs(a(j,i));
%交换第k行和第maxL行begin
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
forj=1:
tmp=a(j,maxC);
a(j,maxC)=a(j,k);
a(j,k)=tmp;
tmp=b_tmp(maxC);
b_tmp(maxC)=b_tmp(k);
b_tmp(k)=tmp;
%交换第k列和第maxC列end
%选主元end
fori=k+1:
n
ifabs(a(k,k))>
1e-6
mult=a(i,k)/a(k,k);
forj=k+1:
n
a(i,j)=a(i,j)-mult*a(k,j);
end%endfor
b(i)=b(i)-mult*b(k);
else
%success=0;
%Gauss消去法失败
disp('
列主元Gauss消去法失败'
);
pause;
exit;
end%endif
end%endfor
%消去过程end
%回代过程begin
x(n)=b(n)/a(n,n);
fori=n-1:
-1:
1
s=0;
forj=i+1:
s=s+a(i,j)*x(j);
x(i)=(b(i)-s)/a(i,i);
end
%回代过程end
%获得正确的解
fori=1:
u(b_tmp(i))=x(i);
End
5.实验总结:
zeros(m,n)生成m行,n列的零矩阵
ones(m,n)生成m行,n列的元素全为1的矩阵
eye(n)生成n阶单位矩阵
rand(m,n)生成m行,n列(0,1)上均匀分布的随机矩阵
diag(x)返回由向量x的元素构成的对角矩阵
tril(A)提取矩阵A的下三角部分生成下三角矩阵
triu(A)提取矩阵A的上三角部分生成上三角矩阵
rank(A)返回矩阵A的秩
det(A)返回方阵A的行列式
inv(A)返回可逆方阵A的逆矩阵
[V,D]=eig(A)返回方阵A的特征值和特征向量
norm(A,p)矩阵或向量的p范数
cond(A,p)矩阵的条件数
[L,U,P]=lu(A)选列主元LU分解
R=chol(X)平方根分解
Hi=hilb(n)生成n阶Hilbert矩阵
通过本次实验基本上掌握并学会用列主元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),
Xi(k+1)=(1/aii)(bi-∑aijxj(k+1)-∑aijxj(k))i=1,2,…,n.
(4)若||x(k+1)-x(k)||2<
则停止计算(x(k+1)作为方程组Ax=b的解)。
(5)若K=N,则停止计算(输出某些信息);
否则置K=K+1,转(3)。
例用Gauss-Seidel迭代法求解下列线性代数方程组,要求||x(k+1)-x(k)||2<
0.0001.
4-10-100x10
-14-10-10x25
0-14-10-1x3=-2
-10-14-10x45
0-10-14-1x5-2
00-10-14x66
四.实验代码:
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'
表示迭代失败
ifnargin<
4max1=100;
3delta=1e-5;
n=length(A);
k=0;
x=zeros(n,1);
y=zeros(n,1);
ok'
;
while1
y=x;
fori=1;
z=b(i);
forj=1:
ifj~=i
z=z-A(i,j)*x(j);
if(abs(A(i,j))<
1e-10)|(k==max1)
flag='
return;
z=z/A(i,i);
x(i)=z;
ifnorm(y-x,inf)<
delta
break;
k=k+1;
五.实验结果:
>
A=[4-10-100;
-14-10-10;
0-14-10-1;
-10-14-10;
0-10-14-1;
00-10-14]
A=
4-10-100
-14-10-10
0-14-10-1
-10-14-10
0-10-14-1
00-10-14
b=[05-25-26]
b=05-25-26
delta=0.0001
delta=1.0000e-004
K>
Gau_Seid(A,b,delta,max1)
11ifnargin<
6.实验总结:
通过这次实验,可分析如下结果:
1)优点分析:
1.程序容错性好,对于能够求解和不能够求解的方程程序都会进行处理并给予回复.
2.程序简单易懂,可读性强。
3.执行效率高,即使是大型矩阵也可以很快求解
2)缺点分析:
1.程序每次只能对一个方程组进行求解,对于矩阵A,b,程序能一次调用就全部求解。
2.算法没能很好的利用MATLAB强大的矩阵运算功能来减少代码,使得代码冗长。
实验四插值法
————三次样条插值法
1.学会用三次样条函数作插值计算
2.通过应用和实际编程计算进一步了解该方法的原理和优缺点。
实验算法(求三次样条插值勤函数S(x)在一列点x1,x2,…xm上的近似值);
(1)输入xi,yi,(i=0,1,…,n),边界条件y0,yn,及xj(j=1,2,…m).
(2)计算hi.hi=xi,xi-1i=1,2,…n.
(3)按下列公式计算μi,λi,gi,i=1,2,…n-1
μi=hi/(hi+hi+1),λi,μi
gi=6/(hi+hi+1)(yi+1-yi)/hi+1-(yi-yi-1)/hi)
(4)按下列公式计算g0,gn,
g0=(6/h1)((y1-y0)/h1-y0`)
gn=6/hn(yn`-(yn-yn-1/hn))
(5)用追赶法求解如下方程组得Mi,i=0,1,…n,
21M0g0
μ12λ1 M1g1
………﹕=﹕
μn-12λn-1Mn-1gn-1
12Mngn
(6)确定点xj(j=1,2,…m)所在的区间[xi-1,xi],并按下列公式计算S(xj);
S(xj)=Mi-1(xi-x)3+Mi((x-xi-1)3/6hi)+(yi-1-(Mi-1/6)hi2)((xi-x)/hi)+(yi-(Mi/6)hi2)((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
处点的纵坐标(要求该曲线具有二阶光滑度)
k
0123456
xk
0.253.18.017.9528.6539.6250.65
yk
5.287949.413..8420.224.928.4431.1
789101112
78104.6156.6208.6260.7312.5
3536.536.634.631.026.34
131415161718
364.4416.3468494507520
20.914.87.83.71.50.2
四.实验代码:
functionS=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);
fork=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);
m(n)=u(n-1)/b(n-1);
fork=n-2:
m(k+1)=(u(k)-c(k)*m(k+2))/b(k);
m
(1)=3*(d
(1)-dx0)/h
(1)-m
(2)/2;
m(n+1)=3*(dxn-d(n))/h(n)-m(n)/2;
fork=0:
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);
五.实验结果
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.17310.20000
-1.01921.24621.27310.5000
0.6558-1.81150.65582.0000
六.实验总结:
plot(x,y)作出以数据(x(i),y(i))为节点的折线图,其中x,y为同长度的向量
subplot(m,n,k)将图形窗口分为m*n个子图,并指向第k幅图
yi=interp1(x,y,xi)根据数据(x,y)给出在xi的分段线性插值结果yi
pp=spline(x,y)返回样条插值的分段多项式(pp)形式结构
pp=csape(x,y,‘边界类型’,‘边界值’)生成各种边界条件的三次样条插值
yi=ppval(pp,xi)pp样条在xi的函数值
ZI=interp2(x,y,z,xi,yi)x,xi为行向量,y,yi为列向量,z为矩阵的双线性二维插值
ZI=interp2(…,'
spline'
)使用二元三次样条插值
ZI=griddata(x,y,z,xi,yi)x,y,z均为向量(不必单调)表示数据,xi,yi为网格向量的三角形线性插值(不规则数据的二维插值)
通过本次实验,考虑一个固定的区间上用插值逼近一个函数。
显然拉格朗日插值中使用的节点越多,插值多项式的次数就越高。
我们自然关心插值多项式的次数增加时,
是否也更加靠近被逼近的函数。
龙格(Runge)给出一个例子是极著名并富有启发性的。
项式插值是不收敛的,即插值的节点多,效果不一定就好。
按一定的规则分别选择等距或者非等距的插值节点,并不断增加插值节点的个数,可以用MATLAB的函数“spline”作此函数的三次样条插值。
实验五数值积分法
1.学会用复合辛普森公式计算某积分的近似值;
2.学会用半步长复合梯形法计算某积分的近似值;
3.通过实际编程计算进一步了解方法的原理和优缺点。
实验算法:
1.复合辛普森公式
∫baf(x)≈h/6[f(a)+4∑f(Xk+1/2)+2∑f(Xk)+f(b)]
其中Xk=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)
其中Xk=a+k*h,(k=0,1,...,n-1),h=(b-a)/n,Xk+1/2=Xk+h/2。
三.实验内容:
1.用复合辛普森公式计算积分∫048(1+cos2x)1/2dx,是误差不超过10-4
2.用半步长复合梯形法的递推算法计算积分
∫01(1-e-x)/xdx
使误差不超过10-4
四.实验代码:
functions=trapr1(f,a,b,n)
h=(b-a)/n;
s=0;
x=a+h*k;
s=s+feval('
lp'
x);
s=h*(feval('
a)+feval('
b))/2+h*s;
functionz=lp(x)
z=1/(1+x^2);
z=(1-exp(-x)^(1/2))/x;
五.实验结果:
trapr1('
0.00000001,1,10)
ans=
0.4439
-1,1,10)
1.5675
diff(x)如果x是向量,返回向量x的差分;
如果x是矩阵,则按各列作差分
diff(x,k)k阶差分
q=polyder(p)求得由向量p表示的多项式导函数的向量表示q
Fx=gradient(F,x)返回向量F表示的一元函数沿x方向的导函数F'
(x),其中x是与F同维数的向量
z=trapz(x,y)x表示积分区间的离散化向量;
y是与x同维数的向量,表示被积函数;
z返回积分的近似值
z=guad(fun,a,b,tol)自适应步长Simpson积分法求得Fun在区间[a,b]上的定积分,Fun为M文件函数句柄,tol为积分精度
z=dblquad(fun,a,b,c,d,tol,method)求得二元函数Fun(x,y)的重积分
z=triplequad(fun,a,b,c,d,e,f,tol,method)求得三元函数Fun(x,y,z)的重积分
通过本次实验可知,线性的积分方程的数值求解,可以被转化为线性代数方程组的求解问题。
而线性代数方程组所含未知数的个数,与用来离散积分的数值方法的节点个数相同。
在节点数相同的前提下,高斯数值积分方法有较高的代数精度,用它通常会得到较好的结果。
高维空间中的积分,如果维数不很高且积分区域是规则的或者能等价地写成
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 报告