2016-2017 第一学期数值分析
上机实验报告
姓名： xxx
学号： 20162…….
学院：土木工程学院
导师：………..
联系电话：…………..
指导老师：………..
目录
第一题 (1)
1.1题目要求 (1)
1.2程序编写 (1)
1.3计算结果及分析 (2)
第二题 (4)
2.1题目要求 (4)
2.2程序编写 (4)
2.3计算结果及分析 (6)
第三题 (7)
3.1题目要求 (7)
3.2程序编写 (7)
3.3计算结果及分析 (8)
第四题 (9)
4.1题目要求 (9)
4.2程序编写 (9)
4.3计算结果及分析 (10)
第五题 (11)
5.1题目要求 (11)
5.2程序编写 (12)
5.3计算结果及分析 (13)
第六题 (17)
6.1题目要求 (17)
6.2程序编写 (17)
6.3计算结果及分析 (18)
6.4程序改进 (18)
第一题
选做的是第（1）小问。
1.1题目要求
编出不动点迭代法求根的程序；把x3+4x2?10=0写成至少四种x=g(x)的形式，取初值x0=1.5,进行不动点迭代求根，并比较收敛性及收敛速度。1.2程序编写
1.3计算结果及分析
①第一种迭代公式：x=x3+4x2+x?10；matlab计算结果如下：（以下为命令行窗口的内容）
2；matlab计算结果如下：
②第二种迭代公式：x=√(10?x3)/4
（以下为命令窗口内容）
2；matlab计算结果如下：
③第三种迭代公式：x=√10
（以下为命令窗口内容）
?；matlab计算结果如下：
④第四种迭代公式：x=10(x2+4x)
（以下为命令窗口内容）
上述4种迭代公式，1、4两种由于在x真实值附近|g`(x)|>1,不满足迭代局部收敛条件，所以迭代序列不收敛。对于2、3两种式子，由于在x真实值附近|g`(x)|<=L<1,满足迭代局部收敛条件,所以迭代序列收敛。对于2、3两迭代公式，由于L3 第二题 选做的是第（2）小问 2.1题目要求 编写有效程序解线性方程组Ax=b。 2.2程序编写 2.3计算结果及分析 （命令行窗口显示内容）数据太多截取一部分 第三题 3.1题目要求 在区间[-1,1]上取x i=?1+0.2i（i=0,1,2?,10)， 对函数f(x)=1 1+25x2 （a）对函数进行多项式插值和三次样条插值，并画出插值函数及f(x)的函数； （b）对函数求其三次拟合曲线并画出拟合曲线的图像，与（a）中结果进行比较。 3.2程序编写 3.3计算结果及分析 matlab画出的图像如下： 从上图可以清晰的看到，三次样条插值曲线与原始曲线最为接近，作为f(x)的插值函数比较理想。而三次多项式拟合的话，与原始结果偏差很大。 第四题 选择的是第一小问 4.1题目要求 编写Gauss -Legendre 求积公式程序，并计算下列积分 （1）∫e ?(cosx)2 dx 1 ?1 （2）∫e ?(cosx)2 dx 2 0 （3）∫x 10|sin?(1x ?)|dx （4）∫ln?(1+x)x ?10dx 4.2程序编写 4.3计算结果及分析 (a) ∫e ?(cosx)2 dx 1?1 （以下为命令行窗口的内容） (b) ∫e ?(cosx)2 dx 2 （以下为命令行窗口的内容） (c) ∫x 1 0|sin? (1x ?)|dx （以下为命令行窗口的内容） (d) ∫log? (1+x)x ?1 0dx （以下为命令行窗口的内容） 第五题5.1题目要求 给定初值问题 f(x)={y`=?50y+50x2+2x,0≤x≤1 y(0)=13???? 用经典的四阶Lunge-Kutta方法求解，步长分别取为?=0.1,0.025,0.01，计算并打印x=0.1i(i=1,2,?,10)各点的值，并与准确解y(x)=1 3 e?50x+x2作比较。 5.2程序编写 5.3计算结果及分析 （以下为命令行窗口的内容）请输入步长h：0.1 从上面的结果可以看到，当步长为0.1时，求解是很不准确的，可能是步长过大导致的，所以应适当减小步长，再尝试进行计算。下面的运算取步长为0.025，再进行试算。 （以下为命令行窗口的内容） 请输入步长h：0.025 从上面的结果可以看到，由于步长h取值减小,误差变得很小，可以看出其误差限是呈线性增加的。 （以下为命令行窗口的内容） 请输入步长h：0.01 从上面的结果可以看到，当步长取0.01时，误差减小了10多倍。此外，我也尝试了采用更小的步长进行计算，计算所得的相对误差随步长减小也越来越小，限于篇幅，也就不再一一罗列计算结果。所以减小步长，可以很好地提高该算法的精确度，若我们把计算结果看成关于步长h的序列，该序列的收敛速度还比较快。而在这一点上，我们也可以通过理论证明，当步长h趋于零时，计算结果就趋于真实结果，即收敛性。 第六题 6.1题目要求 编程计算∑1 cn ∞n=1,其中c =4.4942×10 307 ，给出并观察计算结果，若有问题，分析之。要求逐项相加编程求解，不允许使用matlab 中的现有函数。 6.2程序编写 6.3计算结果及分析 根据我们所学的数学知识，这个关于n的求和序列实际是不收敛的，但是由于计算机计算时的大数吃掉了小数后，数值无法累加。所以显示为上面的收敛结果。 6.4程序改进 上面的算法发生了大数吃小数的现象，为了避免大数吃小数，可以第一步设n个数相加，然后把这些数分为两个一组，每组两个数相加，然后再分为两个一组，循环下去。这样应该能避免大数吃小数。 数值分析实习报告 姓名：gestepoA 学号：201******* 班级：***班 序言 随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。而且还能减少大量的人工计算。 由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JAVA的辅助软件来解决，得到一个满足误差限的解。本文所计算题目，均采用MATLAB进行编程，MATLAB被称为第四代计算机语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。它具有以下优点： 1友好的工作平台和编程环境。MATLAB界面精致，人机交互性强，操作简单。 2简单易用的程序语言。MATLAB是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（M 文件）后再一起运行。 3强大的科学计算机数据处理能力。包含大量计算算法的集合，拥有600多个工程中要用到的数学运算函数。 4出色的图像处理功能，可以方便地输出二维图像，便于我们绘制函数图像。 目录 1 第一题 (4) 1.1 实验目的 (4) 1.2 实验原理和方法 (4) 1.3 实验结果 (5) 1.3.1 最佳平方逼近法 (5) 1.3.2 拉格朗日插值法 (7) 1.3.3 对比 (8) 2 第二题 (9) 2.1实验目的 (9) 2.2 实验原理和方法 (10) 2.3 实验结果 (10) 2.3.1 第一问 (10) 2.3.2 第二问 (11) 2.3.3 第三问 (11) 3 第三题 (12) 3.1实验目的 (12) 3.2 实验原理和方法 (12) 3.3 实验结果 (12) 4 MATLAB程序 (14) 电子科技大学电子工程学院标准实验报告（实验）课程名称MATLAB与数值分析 学生姓名：李培睿 学号：2013020904026 指导教师：程建 一、实验名称 《MATLAB与数值分析》第一次上机实验 二、实验目的 1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算 操作。（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序） 2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号 转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。（用.m文件编写进行符号因式分解和函数求反的程序） 3. 掌握Matlab函数的编写规范。 4、掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、 三维曲线和面的填充、三维等高线等。（用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释） 5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。 三、实验内容 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以x， y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图，允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。不允许使用sort 函数。 四、实验数据及结果分析 题目一： ①在Editor窗口编写函数代码如下： 第2章牛顿插值法实现 参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55. 求牛顿插值多项式和差商的MA TLAB 主程序： function[A,C,L,wcgs,Cw]=newpoly(X,Y) n=length(X);A=zeros(n,n);A(:,1) =Y'; s=0.0;p=1.0;q=1.0;c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1; end C=A(n,n);b=poly(X(n));q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C);Q=poly2sym(q1); syms M wcgs=M*Q/c1;Cw=q1/c1; （1）保存名为newpoly.m 的M 文件 （2）输入MA TLAB 程序 >> X=[242,243,249,250]; >> Y=[13.681,13.526,13.098,13.095]; >> [A,C,L,wcgs,Cw]=newpoly(X,Y) 输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其 ) ()()1(ξ+n n f x R 的系数向量Cw 。 A = 13.6810 0 0 0 13.5260 -0.1550 0 0 13.0980 -0.0713 0.0120 0 13.0950 -0.0030 0.0098 -0.0003 C = 1.0e+003 * 第6章 MATLAB 数值计算 例6.1 求矩阵A 的每行及每列的最大和最小元素，并求整个矩阵的最大和最小元素。 1356 78256323578255631 01-???? -? ?=???? -??A A=[13,-56,78;25,63,-235;78,25,563;1,0,-1]; max(A,[],2) %求每行最大元素 min(A,[],2) %求每行最小元素 max(A) %求每列最大元素 min(A) %求每列最小元素 max(max(A)) %求整个矩阵的最大元素。也可使用命令：max(A(:)) min(min(A)) %求整个矩阵的最小元素。也可使用命令：min(A(:)) 例6.2 求矩阵A 的每行元素的乘积和全部元素的乘积。 A=[1,2,3,4;5,6,7,8;9,10,11,12]; S=prod(A,2) prod(S) %求A 的全部元素的乘积。也可以使用命令prod(A(:)) 例6.3 求向量X =(1！，2！，3！，…，10！)。 X=cumprod(1:10) 例6.4 对二维矩阵x ，从不同维方向求出其标准方差。 x=[4,5,6;1,4,8] %产生一个二维矩阵x y1=std(x,0,1) y2=std(x,1,1) y3=std(x,0,2) y4=std(x,1,2) 例6.5 生成满足正态分布的10000×5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。 X=randn(10000,5); M=mean(X) D=std(X) R=corrcoef(X) 例6.6 对下列矩阵做各种排序。 185412613713-?? ??=?? ??-?? A A=[1,-8,5;4,12,6;13,7,-13]; sort(A) %对A 的每列按升序排序 -sort(-A,2) %对A 的每行按降序排序 [X,I]=sort(A) %对A 按列排序，并将每个元素所在行号送矩阵I 例6.7 给出概率积分 2 (d x x f x x -? e 的数据表如表6.1所示，用不同的插值方法计算f (0.472)。 x=0.46:0.01:0.49; %给出x ，f(x) f=[0.4846555,0.4937542,0.5027498,0.5116683]; format long interp1(x,f,0.472) %用默认方法，即线性插值方法计算f(x) interp1(x,f,0.472,'nearest') %用最近点插值方法计算f(x) interp1(x,f,0.472,'spline') %用3次样条插值方法计算f(x) interp1(x,f,0.472,'cubic') %用3次多项式插值方法计算f(x) format short 例6.8 某检测参数f 随时间t 的采样结果如表6.2，用数据插值法计算t =2，7，12，17，22，17，32，37，42，47，52，57时的f 值。 T=0:5:65; X=2:5:57; MATLAB 实验一 MATLAB 数值计算 试验报告说明： 1 做试验前请先预习，并独立完成试验和试验报告。 2 报告解答方式：将MATLAB 执行命令和最后运行结果从命令窗口拷贝到每题的题目下面，请将报告解答部分的底纹设置为灰色，以便于批阅。 3 在页眉上写清报告名称，学生姓名，学号，专业以及班级。 3 报告以Word 文档书写。 文档命名方式: 学号+姓名+_（下划线）+试验几.doc 如：110400220张三_试验1.doc 4 试验报告doc 文档以附件形式发送到maya_email@https://www.doczj.com/doc/c718043526.html, 。凡文档命名不符合规范，或者发送方式不正确，不予登记。 5 每次试验报告的最后提交期限：下次试验课之前。 一 目的和要求 1 熟练掌握MATLAB 变量的使用 2 熟练掌握矩阵的创建 3 熟练掌握MATLAB 的矩阵和数组的运算 4 使用元胞数组和结构数组 二 试验内容 1 创建矩阵（必做） 1.1使用直接输入，from:step:to ，linspace ，logspace 等方式创建矩阵。 1.2 输入矩阵12342 46836 9 12a ?? ?= ? ?? ? 1.2-1）分别使用全下标和单下标达方式取出元素“8” >>a=[1 2 3 4;2 4 6 8;3 6 9 12] >> a(2,4) %全下标方式 >> a(11) % 单下标方式 1.2-2）分别用不同的方式从矩阵a 中取出子矩阵?? ??? ?1286 4 3 2 %方法一：全下标方式 a([2,3],[1 2 4]) %方法二：单下标方式 a([2 5 11;3 6 12]) % 方法三：利用逻辑向量 l1=logical([0 1 1]) 列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1); end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x') Matlab作业3（数值分析） 机电工程学院（院、系）专业班组 学号姓名实验日期教师评定 1.计算多项式乘法(x2+2x+2)(x2+5x+4)。 答： 2. （1）将(x-6)(x-3)(x-8)展开为系数多项式的形式。（2）求解在x=8时多项 式(x-1)(x-2) (x-3)(x-4)的值。 答：（1） （2） 3. y=sin(x)，x从0到2π，?x=0.02π，求y的最大值、最小值、均值和标准差。 4.设ｘ＝[0.00.30.8 1.1 1.6 2.3]'，y=[0.500.82 1.14 1.25 1.35 1.40]'，试求二次多项式拟合系数，并据此计算x1=[0.9 1.2]时对应的y1。解：x=[0.0 0.3 0.8 1.1 1.6 2.3]'; %输入变量数据x y=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据y p=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数p x1=[0.9 1.2]; %输入点x1 y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318 y1 = a) 1.2909 5.实验数据处理：已知某压力传感器的测试数据如下表 p为压力值，u为电压值，试用多项式 d cp bp ap p u+ + + =2 3 ) ( 来拟 合其特性函数，求出a,b,c,d，并把拟合曲线和各个测试数据点画在同一幅图上。解： >> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9]; u=[10,11,13,14,17,18,22,24,29,34,39]; x=polyfit(p,u,3) %得多项式系数 t=linspace(0,10,100); y=polyval(x,t); %求多项式得值 plot(p,u,'*',t,y,'r') %画拟和曲线 x = 0.0195 -0.0412 1.4469 9.8267 MATLAB 编程题库 1.下面的数据表近似地满足函数2 1cx b ax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像. 625 .0718.0801.0823.0802.0687.0606.0356.0995 .0628.0544.0008.0213.0362.0586.0931.0i i y x ---- 解： x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y; a=z(1); b=z(2); c=z(3); xh=-1:0.1:1; yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*') 2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数 10 的近似根，并写出调用方式： 精度为10 解： >> edit gexianfa.m function [x iter]=gexianfa(f,x0,x1,tol) iter=0; while(norm(x1-x0)>tol) iter=iter+1; x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end >> edit f.m function v=f(x) v=x.*log(x)-1; >> edit g.m function z=g(y) z=y.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8 实验一 误差分析 实验（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（）中19x 的系数作一个小的扰动。我们希望比较（）和（）根的差别，从而分析方程（）的解对扰动的敏感性。 实验内容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数 poly(v)b = 的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve = ))20:1((ve poly roots + 上述简单的Matlab 程序便得到（）的全部根，程序中的“ess ”即是（）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。 如果扰动项的系数ε很小，我们自然感觉（）和（）的解应当相差很小。计算中你有什么出乎意料的发现表明有些解关于如此的扰动敏感性如何 （2）将方程（）中的扰动项改成18x ε或其它形式，实验中又有怎样的现象 出现 （3）（选作部分）请从理论上分析产生这一问题的根源。注意我们可以将 方程（）写成展开的形式， ) 3.1(0 ),(1920=+-= x x x p αα 同时将方程的解x 看成是系数α的函数，考察方程的某个解关于α的扰动是否敏感，与研究它关于α的导数的大小有何关系为什么你发现了什么现象，哪些根关于α的变化更敏感 思考题一：（上述实验的改进） 在上述实验中我们会发现用roots 函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。 1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组， 1231231234748212515 x x x x x x x x x -+=?? -+=-??-++=? （1）给出Jacobi 迭代法的迭代方程，并判定Jacobi 迭代法求解此方程组是否收敛。 （2）若收敛，编程求解该线性方程组。 解（1）：A=[4 -1 1;4 -8 1;-2 1 5] %线性方程组系数矩阵 A = 4 -1 1 4 -8 1 -2 1 5 >> D=diag(diag(A)) D = 4 0 0 0 -8 0 0 0 5 >> L=-tril(A,-1) % A 的下三角矩阵 L = 0 0 0 -4 0 0 2 -1 0 >> U=-triu(A,1) % A 的上三角矩阵 U = 0 1 -1 0 0 -1 0 0 0 B=inv(D)*(L+U) % B 为雅可比迭代矩阵 B = 0 0.2500 -0.2500 0.5000 0 0.1250 0.4000 -0.2000 0 >> r=eigs(B,1) %B 的谱半径 r = 0.3347 < 1 Jacobi迭代法收敛。 (2)在matlab上编写程序如下： A=[4 -1 1;4 -8 1;-2 1 5]; >> b=[7 -21 15]'; >> x0=[0 0 0]'; >> [x,k]=jacobi(A,b,x0,1e-7) x = 2.0000 4.0000 3.0000 k = 17 附jacobi迭代法的matlab程序如下： function [x,k]=jacobi(A,b,x0,eps) % 采用Jacobi迭代法求Ax=b的解 % A为系数矩阵 % b为常数向量 % x0为迭代初始向量 % eps为解的精度控制 max1= 300; %默认最多迭代300，超过300次给出警告D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\(L+U); f=D\b; x=B*x0+f; k=1; %迭代次数 while norm(x-x0)>=eps x0=x; x=B*x0+f; k=k+1; if(k>=max1) disp('迭代超过300次，方程组可能不收敛'); return; end end Tutorial 6 实验报告 实验名称：Matlab数值计算 实验目的： 1、掌握数据统计与分析的方法； 2、掌握数据插值和曲线拟合的方法及其应用； 3、掌握多项式的常用运算。 实验内容： 1.利用randn函数生成符合正态分布的10×5随机矩阵A，进行如下操作： （1）求A的最大元素和最小元素； （2）求A的每行元素的和以及全部元素的和； （3）分别对A的每列元素按升序、每行元素按降序排列。 2.用3次多项式方法插值计算1-100之间整数的平方根。 3.某气象观测站测得某日6：00-18：00之间每隔2h的室内外温度（°C）如下表所示。 使用三次样条插值分别求出该日室内外6：30-17：30之间每隔2h各点的近似温度，并绘制插值后的温度曲线。 4.已知lgx在[1,101]区间10个整数采样点的函数值如下表所示， 试求lgx 的5次拟合多项式p(x)，并绘制lgx 和p(x)在[1,101]区间的函数曲线。 5. 有3个多项式(),(),()P x x x x P x x P x x x =+++=+=++4322 123 24522 3，试进行下列操作： （1） 求()()()()P x P x P x P x =+123。 （2） 求()P x 的根。 （3） 当x 取矩阵A 的每一元素时，求()P x 的值。其中： .....A --?? ??=?? ???? 112140752350525 6. 求函数在指定点的数值导数。 (),,f x x ==123 7. 用数值方法求定积分。 （1）I π =? 210的近似值。 （2）ln() x I dx x += +?1 22011 实验结果： 1. 1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6) disp('牛顿法'); i=1; n0=180; p0=pi/3; tol=10^(-6); for i=1:n0 p=p0-(p0-sin(p0))/(1-cos(p0)); if abs(p-p0)<=10^(-6) disp('用牛顿法求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次牛顿迭代后无法求出方程的解') end 2、disp('Steffensen加速'); p0=pi/3; for i=1:n0 p1=0.5*p0+0.5*cos(p0); p2=0.5*p1+0.5*cos(p1); p=p0-((p1-p0).^2)./(p2-2.*p1+p0); if abs(p-p0)<=10^(-6) disp('用Steffensen加速求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次Steffensen加速后无法求出方程的解') end 1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根， %误差限为 e=10^-4 disp('二分法') a=0.2;b=0.26; tol=0.0001; n0=10; fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1; for i=1:n0 p=(a+b)/2; fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1; if fp==0||(abs((b-a)/2) 实验6 Matlab数值计算 实验目的： 1、掌握数据统计与分析的方法； 2、掌握数据插值和曲线拟合的方法及其应用； 3、掌握多项式的常用运算。 实验内容： 1.利用randn函数生成符合正态分布的10×5随机矩阵A，进行如下操作： （1）求A的最大元素和最小元素； （2）求A的每行元素的和以及全部元素的和； （3）分别对A的每列元素按升序、每行元素按降序排列。 a = randn(10,5)+10; ma = max(max(a)) mi = min(min(a)) s = sum(a,2) sa = sum(sum(a)) p = sort(a) p1 = -sort(-a,2) 2.用3次多项式方法插值计算1-100之间整数的平方根。 f = sqrt(n); interp1(n,f,(1:100),'cubic') 3.某气象观测站测得某日6：00-18：00之间每隔2h的室内外温度（°C）如下表所示。 使用三次样条插值分别求出该日室内外6：30-17：30之间每隔2h 各点的近似温度，并绘制插值后的温度曲线。 n= 6:2:18; f1 = [18 20 22 25 30 28 24]; f2 = [15 19 24 28 34 32 30]; r = 6.5:2:17.5; w = interp1(n,f1,r,'spline'); w1 = interp1(n,f2,r,'spline'); subplot(211),plot(r,w) subplot(212),plot(r,w1) 4. 已知lgx 在[1,101]区间10个整数采样点的函数值如下表所示， 试求lgx 的5次拟合多项式p(x)，并绘制lgx 和p(x)在[1,101]区间的函数曲线。 x = linspace(1,101,10); y = log(x) /log(10); p = polyfit(x,y,5) y1 = polyval(p,x) plot(x,y,':o',x,y1,'-*') legend('sin(x)','fit') 5. 有3个多项式(),(),()P x x x x P x x P x x x =+++=+=++4 3 2 2 123245223，试进 行下列操作： （1） 求()()()()P x P x P x P x =+123。 （2） 求()P x 的根。 （3） 当x 取矩阵A 的每一元素时，求()P x 的值。其中： .....A --?? ? ?=?????? 11214075 2350 5 25 p1 = [1 2 4 0 5]; p2 = [0 0 0 1 2]; roots([1 -1 -1]) x=linspace(0,2*pi,10); y=sin(x); xi=linspace(0,2*pi,100); y1=interp1(x,y,xi); y2=interp1(x,y,xi,'spline'); y3=interp1(x,y,xi,'cublic'); plot(x,y,'o',xi,y1,xi,y2,xi,y3) x=[0 300 600 1000 1500 2000]; y=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; xi=linspace(0,2000,20); yi=1.0332*exp(-(xi+500)/7756); y1=interp1(x,y,xi,'spline'); subplot(2,1,1);plot(x,y,'o',xi,yi,xi,y1,'*') p=polyfit(x,y,2); y2=polyval(p,xi); subplot(2,1,2);plot(x,y,'o',xi,yi,xi,y2,'*') x=[0 300 600 1000 1500 2000]; y=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; xi=linspace(0,2000,20); y1=interp1(x,y,xi,'spline'); subplot(2,1,1);plot(x,y,'-o', xi,y1,'-*') p=polyfit(x,y,2); y2=polyval(p,xi); subplot(2,1,2);plot(x,y,'-o',xi,y2,'-*') 佛山科学技术学院 实 验 报 告 课程名称 数值分析 实验项目 数值积分 专业班级 机械工程 姓 名 余红杰 学 号 2111505010 指导教师 陈剑 成 绩 日 期 月 日 一、实验目的 1、 掌握程序的录入和matlab 的使用和操作； 2、 了解影响线性方程组解的精度的因素——方法与问题的性态。 3、 学会Matlab 提供的“\”的求解线性方程组。 二、实验要求 1、按照题目要求完成实验内容； 2、写出相应的Matlab 程序； 3、给出实验结果(可以用表格展示实验结果)； 4、分析和讨论实验结果并提出可能的优化实验。 5、写出实验报告。 三、实验步骤 1、用LU 分解及列主元高斯消去法解线性方程组 a)??????? ??=??????? ????????? ??----15900001.582012151526099999.2310 7104321x x x x ， 输出b Ax =中系数LU A =分解的矩阵L 和U ，解向量x 和)det(A ;用列主元法的行交换次序解向量x 和求)det(A ；比较两种方法所得结果。 2、用列主高斯消元法解线性方程组b Ax =。 （1）、???? ? ??=????? ??????? ??--11134.981.4987.023.116.427 .199.103.601.3321x x x （2）、???? ? ??=????? ??????? ??--11134.981.4990.023.116.427 .199.103.600.3321x x x 分别输出)det(,,A b A ，解向量x ，（1）中A 的条件数。分析比较（1）、（2）的计算结果 3、线性方程组b Ax =的A 和b 分别为 ??????? ??=1095791068565778710A ,?????? ? ??=31332332b 则解T x ),1,1,1,1(=. 用MATLAB 内部函数求)det(A 和A 的所有特征值和2)(A cond . 若令 ?????? ? ??=+98.99599.6989.998.585604.508.72.71.8710A A δ, 求解b x x A A =++))((δδ，输出向量x δ和2x δ，从理论结果和实际计算两方面分析线性 方程组b Ax =解的相对误差22/x x δ以及A 的相对误差 /A A δ的关系。 四、实验结果 1： %run311.m clc,clear; A = [10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]; b = [8;5.90001;5;1]; %L U 分解 format short %小数点后四位，不然会受到后面的影响 [L U] = lu(A) %解方程组，输出A ，det(A) y = L\b; format long %小数点后15位显示 x = U\y 2015.1.9 上机作业题报告 JONMMX 2000 1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ，其精确值为 )1 1 123(21+--N N 。 （1）编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果 1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序 数值分析幂法与反幂法 matlab程序 随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量。 要求 1）比较不同的原点位移和初值说明收敛性 2）给出迭代结果，生成DOC文件。 3）程序清单，生成M文件。 解答： >> A=rand(5) %随机产生5*5矩阵求随机矩阵 A = 0.7094 0.1626 0.5853 0.6991 0.1493 0.7547 0.1190 0.2238 0.8909 0.2575 0.2760 0.4984 0.7513 0.9593 0.8407 0.6797 0.9597 0.2551 0.5472 0.2543 0.6551 0.3404 0.5060 0.1386 0.8143 >> B=A+A' %A矩阵和A的转置相加，得到随机对称矩阵B B = 1.4187 0.9173 0.8613 1.3788 0.8044 0.9173 0.2380 0.7222 1.8506 0.5979 0.8613 0.7222 1.5025 1.2144 1.3467 1.3788 1.8506 1.2144 1.0944 0.3929 0.8044 0.5979 1.3467 0.3929 1.6286 B=?? ????? ???? ?? ???6286.13929.03467.15979.08044 .03929.00944 .12144.18506 .13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613 .09173 .04187.1 编写幂法、反幂法程序： function [m,u,index,k]=pow(A,u,ep,it_max) % 求矩阵最大特征值的幂法，其中 % A 为矩阵； % ep 为精度要求，缺省为1e-5； % it_max 为最大迭代次数，缺省为100； % m 为绝对值最大的特征值； % u 为对应最大特征值的特征向量； % index ，当index=1时，迭代成功，当index=0时，迭代失败 if nargin<4 it_max=100; end if nargin<3 ep=1e-5; end n=length(A); index=0; k=0; m1=0; m0=0.01; % 修改移位参数，原点移位法加速收敛，为0时，即为幂法 I=eye(n) T=A-m0*I while k<=it_max v=T*u; [vmax,i]=max(abs(v)); m=v(i); u=v/m; if abs(m-m1) 基于MATLAB数值分析实验报告 班级：072115 姓名：李凯 学号：20111003943 实验二：矩阵与向量运算 实验目的：在MATLAB里，会对矩阵与向量进行加、减、数乘、求逆及矩阵特征值运算，以及矩阵的LU分解。 设A是一个n×n方阵，X是一个n维向量，乘积Y=AX可以看作是n维空间变换。如果能够找到一个标量λ，使得存在一个非零向量X，满足：AX=λX （3.1）则可以认为线性变换T（X）=AX将X映射为λX，此时，称X 是对应于特征值λ的特征向量。改写式（3.1）可以得到线性方程组的标准形式：（A-λI）X=0 （3.2）式（3.2）表示矩阵（A-λI）和非零向量X的乘积是零向量，式（3.2）有非零解的充分必要条件是矩阵（A-λI）是奇异的，即：det(A-λI)=0 该行列式可以表示为如下形式： a11–λa12 (1) a21 a22 –λ…a2n =0 (3.3) ………… A n1 a n2 …a nn 将式（3.3）中的行列式展开后，可以得到一个n阶多项式，称为特征多项式： f(λ)=det(A-λI)=(-1)n(λn+c1λn-1+c2λn-2+…+c n-1λ+c n) (3.4) n阶多项式一共有n个根（可以有重根），将每个根λ带入式（3.2），可以得到一个非零解向量。 习题：求下列矩阵的特征多项式的系数和特征值λj： 3 -1 0 A= -1 2 -1 0-1 3 解：在MATLAB中输入命令： A=【3 -1 0；-1 2 -1；0 -1 3】； c=poly(A) roots(c) 得到： 实验四：Lagrange插值多项式 实验目的：理解Lagrange插值多项式的基本概念，熟悉Lagrange插值多项式的公式源代码，并能根据所给条件求出Lagrange插值多项式，理解龙格现象。 %功能：对一组数据做Lagrange插值 %调用格式：yi=Lagran_(x,y,xi) %x,y：数组形式的数据表 %xi：待计算y值的横坐标数组 %yi：用Lagrange还擦之算出y值数组 function fi=Lagran_（x,f,xi） fi=zeros(size(xi)); np1=length(f); for i=1:np1 z=ones(size(xi)); for j=i:np1 if i~=j,z=z.*(xi-x(j))/(x(i)-x(j));end end fi=fi+z*f(i); end return 习题：已知4对数据（1.6，3.3），（2.7，1.22），（3.9，5.61），（5.6，2.94）。写出这四个数据点的Lagrange插值公式，并 Matlab大数值计算题目 1、统计附件1中的数据，对其中的数据划分区间，从0到50，每 10个单位一个区间，分为5个区间，统计每个区间的数量，画出柱状图。 Matlab程序： clear;clc;close all Data=xlsread('数据.xls'); Q=0:10:50; n=length(Data); m=length(Q); T=zeros(size(Q)); for s=1:n for t=1:m-1 if Data(s)>Q(t)&Data(s) 2、统计附件2中第二列数据中1至100每个数字出现的总次数， 附件2中第三列为每出现第二列数字所对应的次数，最后画出柱状图。 Matlab程序： clear;clc;close all Data=load('WEIBOIDWITHCOMMENTS.txt'); DATA=Data(:,2); t=Data(:,3); % m=max(DATA); m=100; T=zeros(m,1); for i=1:m data=DATA; data(data~=ones(size(data))*i)=0; data(data~=0)=1; n=data.*t; N=sum(n); T(i)=N; end bar(T) 3、 找到矩阵迷宫的通路,矩阵第1行第1列为迷宫的入口，第8行 第8列为迷宫的出口。（0表示路，1表示墙） 000 00000011 11010000 01010010 11010010 11010010 00011010 01000011 11110?????????????????????????? Matlab 程序： 主程序： clear all clc maze=[0,0,0,0,0,0,0,0; 0,1,1,1,1,0,1,0; 0,0,0,0,1,0,1,0; 0,1,0,0,0,0,1,0; 0,1,0,1,1,0,1,0; 0,1,0,0,0,0,1,1; 0,1,0,0,1,0,0,0; 0,1,1,1,1,1,1,0]; total=0; maze(1,1)=3;数值分析MATLAB上机实验
《MATLAB与数值分析》第一次上机实验报告
数值分析的matlab实现
第06章_MATLAB数值计算_例题源程序汇总
matlab2012实验1参考答案
数值分析的MATLAB程序
Matlab作业3(数值分析)答案
同济大学数值分析matlab编程题汇编
数值分析实验报告1
matlab数值分析例题
实验6 Matlab数值计算实验报告
数值分析matlab代码
实验6答案 Matlab数值计算
第3章 MATLAB数值计算-习题 答案
MATLAB数值分析实验三(线性方程求解及精度分析)
东南大学-数值分析上机题作业-MATLAB版
数值分析幂法与反幂法-matlab程序
基于MATLAB数值分析实验报告
Matlab大数值计算题目
相关主题
文本预览