数值分析实验讲义.docx
《数值分析实验讲义.docx》由会员分享,可在线阅读,更多相关《数值分析实验讲义.docx(12页珍藏版)》请在冰豆网上搜索。
![数值分析实验讲义.docx](https://file1.bdocx.com/fileroot1/2022-11/23/46c3aefa-27c7-4139-a4b4-481db8778126/46c3aefa-27c7-4139-a4b4-481db87781261.gif)
数值分析实验讲义
《数值分析》实验讲义
卜兵
数值分析实验一
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