数值分析实验讲义.docx

1、数值分析实验讲义数值分析实验讲义卜兵数值分析实验一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 完成以下操作

2、建数组 a=1,2,3,.,20,b=1,3,5,.,39,并求 a, b内积 操作 1) 主窗口点击新建按钮 2) 在弹出的文本编辑窗口添加 a=1:20 b=1:2:39 sum=a*b 3) 单击保存按钮将文件命名为 demo1保存在例1新建文件夹中 4) 在Command Window中输入 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)

3、维 *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=1 2 3;2 3 4;3 4 5; b=1 1 1;2 2 2;3 3 3; a.*b %a,b对应元素相乘 a*b %a,b矩阵相乘 a.b %a对应元素做分母 a./b %b对应元素做分母 * “.”与 b=1 1 1;2 2 2;3 3 3; b3 b.3 b*b*b

4、 %等于b3 例 :编写函数文件demo3实现 sgn函数功能 操作:1)新建 M 文件,并编辑如下 function val=demo3(x) if x0 val=1; elseif xdemo3(0) demo3(90) demo3(-12) 3) 递归调用 例 :编写函数文件demo4,返回输入整数的阶乘 操作:1)新建 M文件,并编辑如下 function val=demo4(n) if n=1|n=0 val=1; else val=n*demo3(n-1); %递归 end 或 function val=demo3(n) val=1; if n=0 val=1; else for

5、i=1:n val=val*i; end end (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

6、,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=

7、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=

8、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=，直到要求的次数位置的各阶展开式在上的图形. 程序代码:function plot_sin_taylor(n)

9、x=-pi:pi/20:pi;plot(x,sin(x);hold on m=fix(n+1)/2); y=0; for k=1:m y=y+(-1)(k-1)/prod(1:2*k-1)*x.(2*k-1); plot(x,y); end axis(-pi pi -1.5 1.5); xlabel(x);ylabel(y);hold off作业:请仿照上例考察与它的各阶Taylor展开式在-1,1内的近似情况. 数值分析实验二（上）第三章 解线性方程组的直接方法1. 实验目的 1)熟悉用Matlab 向量运算 2)用Matlab矩阵求逆及三角分解 3)掌握解线性方程组的直接方法 2. 试验内

10、容 (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); if n(1)=n(2) error(不是方阵!); end n=n(1);det=1; flag=1:n; for k=1:n t=find(abs(A (k:n,k)=max(abs(A (k:n,k);

11、 %寻找主元 t=t(1)+k-1; flag(k)=t; if t=k p=A (k,:); A (k,:)=A (t,:); A (t,:)=p; %交换行 det=-det; end if A (k,k)=0 error(矩阵不可逆!); end det=det*A (k,k); h=1/A (k,k); A (k,k)=h; A (1:k-1 k+1:n,k)=A (1:k-1 k+1:n,k)*(-h); for i=1:n if i=k A (i,1:k-1 k+1:n)=A (i,1:k-1 k+1:n)+A (k,1:k-1. k+1:n)*A (i,k); %消去 end e

12、nd A (k,1:k-1 k+1:n)=A (k,1:k-1 k+1:n)*h; end for k=n-1:-1:1 t=flag(k); if k=t %调整 A 列(因为原矩阵换行) p=A (:,t); A (:,t)=A (:,k); A (:,k)=p; end end 例1: 解方程组 并验证. 解: A=0.6428 0.3475 -0.8468;0.3475 1.8423 0.4759;-0.8468 0.4759 1.2147; b=0.4127;1.7321;-0.8621; x=GaJo_inv(A)*b A*x-b (2) Gauss 消去法变形-选主元的三角变形法

13、 通过交换矩阵的行实现矩阵的 PA=LU 分解,采用与列主元相似的方法, 避免一些奇异情况的情况下分解无法进行.具体实现如下: function L,A,P=LU_P(A) %选主元的三角变形法P181 %L,A,P=LU_P(A) % A 待分解阵 % L 返回单位下三角阵 % A 返回上三角阵,存储在原矩阵中 % P 返回排列阵 n=size(A); if n(1)=n(2) error(不是方阵!); end n=n(1);flag=1:n; P=eye(n,n);L=P; for k=1:n if k=1 A (k:n,k)=A (k:n,k)-A (k:n,1:k-1)*A (1:k

14、-1,k); end t=find(abs(A (k:n,k)=max(abs(A (k:n,k); t=t(1)+k-1; %选主元 flag(k)=t; if t=k p=A (k,:); A (k,:)=A (t,:); A (t,:)=p; end %交换行 if abs(A (k,k) 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_S

15、olve(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; if nargin4 dalt=1e-8; end if nargin=dalt Y=b-A*x; e=max(abs(Y-x); x=Y; i=i+1; end if nargout1 n=i; end 例1 对不同初值用 Jacobi 迭代法解例 1

16、 中方程组AX=b,比较结果. x0=1.5 2.0 3 1 1.5 1.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 初始值

17、,默认为0 % dalt 精度,默认为 10 % x 返回跌代结果 % n 返回跌代次数 r=size(b); a=b; if nargin5 dalt=1e-8; end if nargindalt e=0; for i=1:r t=x(i); x(i)=(1-w)*x(i)+w*(b(i)-A (i,:)*x); t=abs(x(i)-t); if te e=t; end end m=m+1; end if nargout1 n=m; end 例2 对不同松弛因子用SOR解例 1 中方程组A X =b,并与例 2 比较结果. 解: x,n=SOR_Solve(A,f,1.42) x,n=S

18、OR_Solve(A,f,1) 数值分析实验三第四章 函数的数值逼近一. 实验目的 1) Matlab中多项式的表示及多项式运算 2) 用Matlab实现拉格朗日及牛顿插值法 3) 用多项式插值法拟合数据 二. 实验内容 1. 多项式运算 (1) 多项式表示对多项式，用以下的行向量表示：。这样把多项式转化为向量运算。 p=1,-5,6,-33; poly2sym(p) %符号表示多项式 ans = x3-5*x2+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=

19、-5 -3+4i -3-4i; p=poly(root) %生成以 root中元素为根的多项式 p = 1 11 55 125 poly2sym(p) ans =x3+11*x2+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 = 192 416 800 1392 2240 polyval(p,b) %矩阵b对应元素多项式值 ans = 192 287 287 192 polyvalm(p,b) %按矩阵多项式计算 ans = 24

20、8 168 168 248 roots(p) %多项式p的所有根 ans = -5.0000 -3.0000 + 4.0000i -3.0000 - 4.0000i p=2 -5 6 -1 9; d=3 -90 -18; pd=conv(p,d) %多项式p,d相乘 pd = 6 -195 432 -453 9 -792 -162 p1=deconv(pd,d) %p1 为多项式pd 除 d(等于p) p1 = 2 -5 6 -1 9 p=1 2 3 4 5 6; q=3 2 1; s,r=deconv(p,q) %s 为商,r余项 s = 0.3333 0.4444 0.5926 0.790

21、1 r = 0 0 0 0 2.8272 5.2099 Dp=polyder(p) %对p求导所得多项式 Dp = 5 8 9 8 5 2拉格朗日插值法 (1) 算法实现 function p=Lag_polyfit(X,Y) %Matlab函数文件 Lag_polyfit.m %p=Lag_polyfit(X,Y) % X 拟合自变量 % Y 拟合函数值 % p 所得的拟合多项式系数 if size(X)=size(Y) error(变量不匹配); end %如果要拟合的函数值与自变量维数不一样,则退出报错 tic %开始记时 format long g %设置最合适的数字格式 r=size

22、(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=1 1; r =A ; %A,r 为长度为 2的(1,2)向量 for i=1:n %计算循环开始 A=1,-X (i); %A 为一次多项式x-x(i)系数 p0,r=deconv(W,A); %进行多项式除法W/A,p0 为商 b=Y (i)/z(i); p0

23、=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 函数(蓝色)与多项式(

24、红色)效果图 命令窗口中输入 lagsin观察图形输出窗口 lagsin %函数(文件)名不区分大小写 3牛顿插值法 (1) 均差表计算及算法实现 function p=Newton_Polyfit(X,Y) %Matlab函数文件 Newton_polyfit.m %牛顿插值法多项式拟合 %p=Newton_Polyfit(X,Y) % X 拟合自变量 % Y 拟合函数值 % p 所得的拟合多项式系数 if size(X)=size(Y) error(变量不匹配?);%如果要拟合的函数值与自变量维数不一样,则退出报错 end format long g %设置最合适的数字格式 r=size(

25、X); n=r(2); %n 为要拟合的数据长度 M=ones(n,n); %均差表工作矩阵 M(:,1)=Y; %设置均差表第一列为函数值 for i=2:n for j=i:n %高阶均差推算 M(j,i)=(M(j,i-1)-M(j-1,i-1)/ (X (j)-X (j-i+1); end end M %显示均差表 p0=zeros(1,n-1) M(1,1);p=p0 ; %拟合多项式存储向量及初值 for i=1:n-1 p1=M(i+1,i+1).*poly(X (1:i); %对应阶均差的多项式 p0=zeros(1,n-i-1) p1; p=p+p0 ; %对应项累加 end

26、 (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 数值分析实验四第五章

27、 数值积分 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); if m=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