数值分析实验讲义.docx
《数值分析实验讲义.docx》由会员分享，可在线阅读，更多相关《数值分析实验讲义.docx（12页珍藏版）》请在冰豆网上搜索。
数值分析实验讲义
《数值分析》实验讲义
卜兵
数值分析实验一
Matlab基本操作
1.学习目的
1）熟悉Matlab的运行环境及各种窗口2）掌握Matlab的矩阵变量类型,矩阵输入和矩阵的基本运算3）掌握命令及函数文件的作用及区别,并编写简单的M文件4）能熟练的向查寻目录中添加新目录,掌握常用的Matlab系统命令
5）掌握Matlab的控制语句
6）熟悉数组运算
7）Matlab图形处理功能
8）Matlab程序初步设计
2.学习内容
（1）Matlab启动与环境设置
（2）Matlab基本运算操作（3）Matlab的文件例1:
编写命令文件demo1完成以下操作建数组a=[1,2,3,...,20],b=[1,3,5,...,39],并求a,b内积操作1）主窗口点击新建按钮
2）在弹出的文本编辑窗口添加a=1:
20b=1:
2:
39sum=a*b'3）单击保存按钮
将文件命名为demo1保存在例1新建文件夹中4）在CommandWindow中输入demo1并回车
（4）数组运算（相同类型的运算）
1）’:
’引用*A（:
n）矩阵A的n列所有元素>>A=rand（4,5）;>>A（:
3）=（1:
4）’%引用的为一列向量*A（m,:
）矩阵A的m行所有元素>>A（4,:
）=2:
6*A（:
）矩阵A所有元素
>>A（:
）2）维*reshape（X,M,N,P,..）将已知矩阵X为M*N*P..矩阵>>a=1:
12;>>b=reshape（a,2,6）*用’:
’引用>>a=zeros（3,4）;>>a（:
）=1:
12%Matlab矩阵元素按列存储>>a（4）>>a（1,2）3）‘.’运算同类型矩阵元素对应元素运算*“.*”,“./”与”.\’运算>>a=[123;234;345];>>b=[111;222;333];>>a.*b%a,b对应元素相乘>>a*b%a,b矩阵相乘>>a.\b%a对应元素做分母>>a./b%b对应元素做分母*“.^”与^>>b=[111;222;333];>>b^3>>b.^3>>b*b*b%等于b^3
例:
编写函数文件demo3实现sgn函数功能操作:
1）新建M文件,并编辑如下functionval=demo3（x）ifx>0val=1;elseifx<0val=-1;elseval=0;end2）将文件保存在查询目录内3）>>demo3（0）>>demo3（90）>>demo3（-12）3）递归调用例:
编写函数文件demo4,返回输入整数的阶乘操作:
1）新建M文件,并编辑如下functionval=demo4（n）ifn==1|n==0val=1;elseval=n*demo3（n-1）;%递归end或function[val]=demo3（n）val=1;ifn==0val=1;elsefori=1:
nval=val*i;endend
（6）Matlab图形处理初步
1）二维图形plot（x,y,s）例1:
>>x=rand（100,1）;>>y=rand（100,1）;>>z=x+y.*i;>>plot（y）>>plot（z）;例2:
>>x=0.1:
0.01*pi:
pi;>>y=sin（x）.*cos（x）;>>plot（x,y）;注意:
当两个输入变量同为向量时,x,y维数相同.x,y为同阶矩阵时将按列或行进行.例3:
>>x=0.1:
0.01*pi:
pi;>>y=[sin（x）',cos（x）'];>>plot（[x'],y）>>plot（[x',x'],y）>>plot（x',y（:
1）,x',y（:
2））例4:
>>x=0.1:
0.1*pi:
2*pi;>>y=sin（x）;>>z=cos（x）;>>plot（x,y,'--k',x,z,'-.rd'）注:
s图形设置选项选项说明选项说明-实线y黄色:
点线r红色-.点划线g绿色..虚线k黑色o圆号++号**号d菱形
2）三维图形*plot3（x,y,z,s）%其中x,y和z为3个相同维数的向量*plot3（X,Y,Z,s）%其中X,Y和Z为3个相同阶数的矩阵,函
数绘3矩阵的列向量曲线*plot3（x1,y1,z1,s1,x2,y2,z2,s2,…）例1:
>>x=0:
pi/50:
10*pi;>>y=sin（x）;>>z=cos（x）;>>plot3（x,y,z）;例2:
>>[x,y]=meshgrid（-2:
0.1:
2,-2:
0.1:
2）;%产生网格点>>z=x.*exp（-x.^2-y.^2）;>>plot3（x,y,z）;例3:
>>x=-8:
0.5:
8;y=x';>>a=ones（size（y））*x;>>b=y*ones（size（x））;>>c=sqrt（a.^2+b.^2）+eps;>>z=sin（c）./c;>>surf（x,y,z）>>mesh（x,y,z）
3）特殊的三维图形函数*[x,y,z]=sphere（n）%画球,n默认值20例1:
>>[a,b,c]=sphere（40）;>>surf（a,b,c）>>axis（'equal'）;>>axis（'square'）;*[x,y,z]=cylinder（R,N）%R母线向量,N分格数例2:
>>x=0:
pi/20:
pi*3;>>r=5+cos（x）;>>[a,b,c]=cylinder（r,30）;>>mesh（a,b,c）
练习：
例：
需要考察函数sin（x）与它的各阶Taylor展开式
的近似情况。
请编写一个函数（程序）。
当输入Taylor展开式的阶数后，函数将自动画出y=sin（x），y=x,y=
，直到要求的次数位置的各阶展开式在[
]上的图形.
程序代码:
functionplot_sin_taylor（n）
x=-pi:
pi/20:
pi;
plot（x,sin（x））;
holdon
m=fix（（n+1）/2）;
y=0;
fork=1:
m
y=y+（-1）^（k-1）/prod（[1:
2*k-1]）*x.^（2*k-1）;
plot（x,y）;
end
axis（[-pipi-1.51.5]）;
xlabel（'x'）;ylabel（'y'）;
holdoff
作业:
请仿照上例考察
与它的各阶Taylor展开式在[-1,1]内的近似情况.
数值分析实验二（上）
第三章解线性方程组的直接方法
1.实验目的1）熟悉用Matlab向量运算2）用Matlab矩阵求逆及三角分解3）掌握解线性方程组的直接方法2.试验内容
（1）Gauss_Jordan列主元消去法求方阵逆.
选择列主元,消去对角线上方和下方的元素.用其实现矩阵求逆如下:
function[A,det]=GaJo_inv（A）%Gauss-Jordan列主元消去法求方阵逆P178%[A,det]=GaJo_inv（A）%A要求逆矩阵%det按需求返回A的行列式%A返回逆矩阵,放在A中n=size（A）;ifn
（1）~=n
（2）error（'不是方阵!
'）;endn=n
（1）;det=1;flag=1:
n;fork=1:
nt=find（abs（A（k:
n,k））==max（abs（A（k:
n,k））））;%寻找主元t=t
（1）+k-1;flag（k）=t;ift~=kp=A（k,:
）;A（k,:
）=A（t,:
）;A（t,:
）=p;%交换行det=-det;endifA（k,k）==0error（'矩阵不可逆!
'）;enddet=det*A（k,k）;h=1/A（k,k）;A（k,k）=h;A（[1:
k-1k+1:
n],k）=A（[1:
k-1k+1:
n],k）*（-h）;fori=1:
nifi~=kA（i,[1:
k-1k+1:
n]）=A（i,[1:
k-1k+1:
n]）+A（k,[1:
k-1….k+1:
n]）*A（i,k）;%消去endendA（k,[1:
k-1k+1:
n]）=A（k,[1:
k-1k+1:
n]）*h;endfork=n-1:
-1:
1t=flag（k）;ifk~=t%调整A列（因为原矩阵换行）p=A（:
t）;A（:
t）=A（:
k）;A（:
k）=p;endend
例1:
解方程组
并验证.
解:
>>A=[0.64280.3475-0.8468;0.34751.84230.4759;-0.84680.47591.2147];>>b=[0.4127;1.7321;-0.8621];>>x=GaJo_inv（A）*b>>A*x-b
（2）Gauss消去法变形----选主元的三角变形法通过交换矩阵的行实现矩阵的PA=LU分解,采用与列主元相似的方法,避免一些奇异情况的情况下分解无法进行.具体实现如下:
function[L,A,P]=LU_P（A）%选主元的三角变形法P181%[L,A,P]=LU_P（A）%A待分解阵%L返回单位下三角阵%A返回上三角阵,存储在原矩阵中%P返回排列阵n=size（A）;ifn
（1）~=n
（2）error（'不是方阵!
'）;endn=n
（1）;flag=1:
n;P=eye（n,n）;L=P;fork=1:
nifk~=1A（k:
n,k）=A（k:
n,k）-A（k:
n,1:
k-1）*A（1:
k-1,k）;endt=find（abs（A（k:
n,k））==max（abs（A（k:
n,k））））;t=t
（1）+k-1;%选主元flag（k）=t;ift~=kp=A（k,:
）;A（k,:
）=A（t,:
）;A（t,:
）=p;end%交换行ifabs（A（k,k））'）;endifk~=nA（k+1:
n,k）=（1/A（k,k））*A（k+1:
n,k）;A（k+1:
n,k）=A（k+1:
n,k）;A（k,k+1:
n）=A（k,k+1:
n）-A（k,1:
k-1）*A（1:
k-1,k+1:
n）;endendfork=1:
n-1L（k+1:
n,k）=A（k+1:
n,k）;A（k+1:
n,k）=zeros（n-k,1）;end%得到L,Ufork=n-1:
-1:
1t=flag（k）;ifk~=t%按主元配制排列阵p=P（:
t）;P（:
t）=P（:
k）;P（:
k）=p;endend例2:
产生4,4随机矩阵,做PLU分解,并作验证.
>>A=rand（4,4）*10;>>[L,U,P]=LU_P（A）>>P*A-L*U
数值分析实验二（下）
第三章解线性方程组的迭代方法
1.实验目的
用Matlab实现Jacobi及超松弛跌代法
2.实验内容
（1）Jacobi迭代法解线性方程组迭代公式:
function[x,n]=Jacobi_Solve（A,b,x0,dalt）%Jacobi跌代法解线性方程组P213%[x,n]=Jacobi_Solve（A,b,x0,dalt）%A方程组系数%b常数项（列向量）%x0初始值,默认为0%dalt精度,默认为10%x返回跌代结果%n返回跌代次数e=1;i=0;r=size（b）;a=b;ifnargin<4dalt=1e-8;endifnargin<3x=zeros（r）;elsex=x0;endr=r
（1）;fort=1:
ra（t）=A（t,t）;A（t,t）=0;A（t,:
）=A（t,:
）/a（t）;endb=b./a;whilee>=daltY=b-A*x;e=max（abs（Y-x））;x=Y;i=i+1;endifnargout>1n=i;end
例1对不同初值用Jacobi迭代法解例1中方程组AX=b,比较结果.>>x0=[1.52.0311.51.8]';>>[x,n]=Jacobi_Solve（A,f',x0）>>[x,n]=Jacobi_Solve（A,f'）
（2）超松弛跌代法解线性方程组
超松弛跌代法是Gauss-Seidel跌代法的一种加速,是解大型稀疏矩阵方程组的有效方法之一.但需要选择好的加速因子（最佳松弛因子）.
function[x,n]=SOR_Solve（A,b,w,x0,dalt）%超松弛跌代法解线性方程组P213%[x,n]=SOR_Solve（A,b,w,x0,dalt）%A方程组系数%b常数项（列向量）%w松弛因子%x0初始值,默认为0%dalt精度,默认为10%x返回跌代结果%n返回跌代次数r=size（b）;a=b;ifnargin<5dalt=1e-8;endifnargin<4x=zeros（r）;elsex=x0;endr=r
（1）;m=0;e=1;fort=1:
ra（t）=A（t,t）;A（t,t）=0;A（t,:
）=A（t,:
）/a（t）;endb=b./a;whilee>dalte=0;fori=1:
rt=x（i）;x（i）=（1-w）*x（i）+w*（b（i）-A（i,:
）*x）;t=abs（x（i）-t）;ift>ee=t;endendm=m+1;endifnargout>1n=m;end
例2对不同松弛因子用SOR解例1中方程组AX=b,并与例2比较结果.
解:
>>[x,n]=SOR_Solve（A,f',1.42）
>>[x,n]=SOR_Solve（A,f',1）
数值分析实验三
第四章函数的数值逼近
一.实验目的
1）Matlab中多项式的表示及多项式运算2）用Matlab实现拉格朗日及牛顿插值法3）用多项式插值法拟合数据
二.实验内容
1.多项式运算
（1）多项式表示
对多项式
，用以下的行向量表示：
。
这样把多项式转化为向量运算。
①>>p=[1,-5,6,-33];>>poly2sym（p）%符号表示多项式ans=x^3-5*x^2+6*x-33②>>a=[1,2,3;2,3,4;3,4,5];>>p1=poly（a）%矩阵a的特征多项式p1=1.0000-9.0000-6.0000-0.0000③>>root=[-5-3+4i-3-4i];>>p=poly（root）%生成以root中元素为根的多项式p=11155125>>poly2sym（p）ans=
x^3+11*x^2+55*x+125
（2）多项式运算①>>p=[1,11,55,125];%建立多项式p>>a=1:
2:
10;%向量a>>b=[1,2;2,1];%方阵b>>polyval（p,a）%计算a对应元素的多项式值ans=19241680013922240>>polyval（p,b）%矩阵b对应元素多项式值ans=192287287192
>>polyvalm（p,b）%按矩阵多项式计算
ans=
248168168248>>roots（p）%多项式p的所有根
ans=
-5.0000
-3.0000+4.0000i
-3.0000-4.0000i
②>>p=[2-56-19];
>>d=[3-90-18];
>>pd=conv（p,d）%多项式p,d相乘
pd=
6-195432-4539-792-162
>>p1=deconv（pd,d）%p1为多项式pd除d（等于p）
p1=
2-56-19
>>p=[123456];
>>q=[321];
>>[s,r]=deconv（p,q）%s为商,r余项
s=
0.33330.44440.59260.7901
r=
00002.82725.2099
>>Dp=polyder（p）%对p求导所得多项式
Dp=
58985
2．拉格朗日插值法
（1）算法实现
function[p]=Lag_polyfit（X,Y）%Matlab函数文件Lag_polyfit.m%[p]=Lag_polyfit（X,Y）%X拟合自变量%Y拟合函数值%p所得的拟合多项式系数ifsize（X）~=size（Y）error（'变量不匹配'）;end%如果要拟合的函数值与自变量维数不一样,则退出报错tic%开始记时formatlongg%设置最合适的数字格式r=size（Y）;n=r
（2）;%n为要拟合的数据长度p=zeros（1,n）;%保存所得多项式系数p0=p;b=0;%工作变量W=poly（X）;%W为以X为根的多项式dW=polyder（W）;%dW为对多项式W求导后的多项式系数z=polyval（dW,X）;%z为以dW为系数的多项式对X的值A=[11];r=A;%A,r为长度为2的（1,2）向量fori=1:
n%计算循环开始A=[1,-X（i）];%A为一次多项式x-x（i）系数[p0,r]=deconv（W,A）;%进行多项式除法W/A,p0为商b=Y（i）/z（i）;p0=b.*p0;%p0为累加项p=p+p0;end%循环结束
disp（toc）%显示所用时间
（2）例题%建立Matlab命令文件LagSin.m%用拉格朗日插值法拟合[0,2*pi]上的sin函数程序命令行：
x0=linspace（0,2*pi,10）;y0=sin（x0）;%拟合数据采样p=Lag_polyfit（x0,y0）;%调用Lag_polyfit计算拟合多项式x=0:
0.2:
2*pi;
y1=sin（x）;%sin函数值y2=polyval（p,x）;%拟合多项式值plot（x,y1,x,y2,'r'） %sin函数（蓝色）与多项式（红色）效果图命令窗口中输入lagsin观察图形输出窗口>>lagsin%函数（文件）名不区分大小写
3．牛顿插值法
（1）均差表计算及算法实现function[p]=Newton_Polyfit（X,Y）
%Matlab函数文件Newton_polyfit.m%牛顿插值法多项式拟合%[p]=Newton_Polyfit（X,Y）%X拟合自变量%Y拟合函数值%p所得的拟合多项式系数ifsize（X）~=size（Y）error（'变量不匹配?
’）;%如果要拟合的函数值与自变量维数不一样,则退出报错endformatlongg%设置最合适的数字格式r=size（X）;n=r
（2）;%n为要拟合的数据长度M=ones（n,n）;%均差表工作矩阵M（:
1）=Y';%设置均差表第一列为函数值fori=2:
nforj=i:
n%高阶均差推算M（j,i）=（M（j,i-1）-M（j-1,i-1））/（X（j）-X（j-i+1））;endendM%显示均差表p0=[zeros（1,n-1）M（1,1）];p=p0;%拟合多项式存储向量及初值fori=1:
n-1p1=M（i+1,i+1）.*poly（X（1:
i））;%对应阶均差的多项式p0=[zeros（1,n-i-1）p1];p=p+p0;%对应项累加end
（2）例题%Matlab命令文件NewtonSin.m%用牛顿插值法拟合[0,2*pi]上的sin函数x0=linspace（0,2*pi,10）;y0=sin（x0）;%拟合数据采样p=Newton_polyfit（x0,y0）;%调用Newton_polyfit计算拟合多项式x=0:
0.2:
2*pi;y1=sin（x）;%sin函数值y2=polyval（p,x）;%拟合多项式值plot（x,y1,x,y2,'r'）%sin函数（蓝色）与多项式（红色）效果图命令窗口中输入newtonsin观察图形输出窗口>>newtonsin
数值分析实验四第五章数值积分
1.实验目的1）熟悉Matlab矩阵操作2）用Matlab实现数值积分Cores,Romberg及复化Gauss算法3）学会数值积分有关应用2.实验内容
（1）Cores积分在使用牛顿-柯斯特公式时,通过提高阶为提高精度的途径并不能取得满意的效果.通常采用复化积分法.以复化柯斯特公式为例:
将每个区间4等分,依此用柯斯特公式function[val]=Cores_int（Y,u）%[val]=Cores_int（Y,u）%Y积分函数数值%u区间间隔%val返回积分值r=size（Y）;n=r
（2）;v=n;h=4*u;m=mod（n-1,4）;ifm==0%如果取值点恰巧满足复化区间数t=（n-1）/4;x=ones（t,1）;val=（h/90）*（7*Y
（1）+32*（Y（2:
4:
n-3）*x）+12*（Y（3:
4:
n-2）*x）+32*（Y（4:
4:
n-1）*x）+….
+14*（Y（5:
4:
n-4）*ones（size（（5:
4:
n-4）'）））+7*Y（n））;e
升级会员 