计算方法实验指导书1.docx
- 文档编号:25617557
- 上传时间:2023-06-10
- 格式:DOCX
- 页数:36
- 大小:212.98KB
计算方法实验指导书1.docx
《计算方法实验指导书1.docx》由会员分享,可在线阅读,更多相关《计算方法实验指导书1.docx(36页珍藏版)》请在冰豆网上搜索。
计算方法实验指导书1
第一章绪论
一、主要要求
通过实验,认真理解和体会数值计算的稳定性、精确性与步长的关系。
二、主要结果回顾:
1、算法:
电子计算机实质上只会做加、减、乘、除等算术运算和一些逻辑运算,由这些基本运算及运算顺序规定构成的解题步骤,称为算法.它可以用框图、算法语言、数学语言或自然语言来描述。
用计算机算法语言描述的算法称为计算机程序。
(如c—语言程序,c++语言程序,Matlab语言程序等)。
2、最有效的算法:
应该运算量少,应用范围广,需用存储单元少,逻辑结构简单,便于编写计算机程序,而且计算结果可靠。
3、算法的稳定性:
一个算法如果输入数据有误差,而在计算过程中舍入误差不增长,则称此算法是数值稳定的,否则称此算法为不稳定的。
换句话说:
若误差传播是可控制的,则称此算法是数值稳定的,否则称此算法为不稳定的。
4、控制误差传播的几个原则:
1)防止相近的两数相减;
2)防止大数吃小数;
3)防止接近零的数做除数;
4)要控制舍入误差的累积和传播;
5)简化计算步骤,减小运算次数,避免误差积累。
三、数值计算实验(以下实验都需利用Matlab软件来完成)
实验1.1(体会数值计算精度与步长关系的实验)
实验目的:
数值计算中误差是不可避免的,要求通过本实验初步认识数值分析中两个重要概念:
截断误差和舍入误差,并认真体会误差对计算结果的影响。
问题提出:
设一元函数f:
R→R,则f在x0的导数定义为:
实验内容:
根据不同的步长可设计两种算法,计算f在x0处的导数。
计算一阶导数的算法有两种:
(1)
(2)
请给出几个计算高阶导数的近似算法,并完成如下工作:
1、对同样的h,比较
(1)式和
(2)式的计算结果;
2、针对计算高阶导数的算法,比较h取不同值时
(1)式和
(2)式的计算结果。
实验要求:
选择有代表性的函数f(x)(最好多选择几个),利用Matlab提供的绘图工具画出该函数在某区间的导数曲线f(s)(x),再将数值计算的结果用Matlab画出来,认真思考实验的结果。
实验分析:
不论采用怎样的算法,计算结果通常都会有误差。
比如算法
(1),由Taylor公式,知:
所以有
利用上式来计算f’(x0),其截断误差为:
所以误差是存在的,并且当步长h越来越小时,
(1)的近似程度也越来越好。
类似地可以分析
(2)的截断误差为:
上述截断误差的分析表明
(2)是比
(1)更好的算法,因为对步长h(<<1),
(2)比
(1)更接近于f’(x0)。
计算方法的截断误差是数值计算中误差的重要来源,但不是唯一的!
!
如果在实验中确已将h取到足够小的话,特别在高阶导数的计算中,就会发现当h小到一定程度后,计算结果的误差不但不再减少,反而会变大!
(参见图1)
事实上当步长h过小时,计算结果的误差变大就是由于舍入误差的缘故。
截断误差用我们原有的数学思维方式就比较容易理解的,而舍入误差则是本课程引入的一个新概念。
要真正理解舍入误差,特别是它在计算中的传播及最终对计算结果的影响,是初步具备科学计算能力的重要标志。
希望大家在完成实验后,认真仔细去体会截断误差和舍入误差的含义及对计算结果的影响。
0
10-7
10-5
10-5
10-4
10-3
10-2
10-1
10-3
10-1
(图1)
实验1.2(理解误差传播与算法稳定性实验)
实验目的:
体会算法的稳定性在选择算法中的地位。
问题提出:
考虑一个简单的积分序列
n=1,2,….
显然In>0,n=1,2,…当n=1时,得:
当n≥2时,由分部积分可得:
n=2,3,……
另外,还有:
实验内容:
由递推关系In=1-nIn-1,可得计算积分序列{In}的两种算法:
算法一、直接使用递推公式得:
In=1-nIn-1n=2,3…
算法二、利用递推公式变形得:
实验要求:
用上述两种算法分别在计算中采用5位、6位和7位有效数字,请判断哪种算法给出的结果更精确。
实验分析:
两种算法的优劣可能与你的第一感觉完全不同。
设算法一中I1的误差为e1,由I1递推计算In的误差为en
算法二中IN的误差为εN,由IN向前递推计算In(n 如果假定上述两种算法在后面的计算中都不再引入其他误差,则有: 由此可见,算法一中的e1可能很小,但在计算中它的影响急剧扩散,结果真实的数据很快被淹没了;而算法二中的εN可能相对比较大,但在计算中误差影响不扩散,某一步产生误差后,该误差对后面的影响不断衰减。 注: 误差扩散的算法是不稳定的,是我们所不期望的;误差衰减的算法是稳定的,是我们努力寻求的,也是贯穿本课程始终的目标。 第二章一元非线性方程的解法 一、主要要求 编写二分法、Newton迭代法和快速弦截法通用子程序。 二、主要结果回顾 1、二分法的基本思想: 二分法就是将方程的有根区间对分,然后再选择比原区间缩小一半的有根区间,如此继续下去,直到得到满足精度要求的根为止的一种简单的区间方法。 基本法原理: 给定方程f(x)=0,设f(x)在区间[a,b]连续,且f(a)×f(b)<0,则方程f(x)在(a,b)内至少有一根,为便于讨论,不妨设方程f(x)=0在(a,b)内只有一实根x*采取使有根区间逐步缩小,从而得到满足精度要求的实根x*的近似值xk。 取[a,b]区间二等分的中点 ,若f(x0)=0,则x0是f(x)=0的实根; 若f(a)f(x0)<0成立,则x*必在区间(a,x0)内,取a1=a,b1=x0;否则x*必在区间(x0,b)内,取a1=x0,b1=b,这样,得到新区间(a1,b1),其长度为[a,b]的一半,如此继续下去,进行k次等分后,得到一系列有根区间: ,其中每一个区间长度都是前一个区间长度的一半,从而[ak,bk]的长度为 如此继续下去,则有这些区间将收敛于一点,该点即为所求的根.当做到第k步时,有: (ε为给定的精度) 此时 即为所求方程的近似解。 2、二分法算法: 1)输入有根区间的端点a,b及预先给定的精度ε; 2)计算(a+b)/2并将其赋值给变量c,记为(a+b)/2→c; 3)若f(a)f(c)<0,则c→b,转向4),否则c→a,转向4); 4)若b-a<ε,则输出满足精度要求的跟c,否则转向2) 3、二分法程序框图(见图2) 4、不动点迭代法: 不动点迭代法又称简单迭代法,基本思想是构造不动点方程,以求得近似根。 即由方程f(x)=0变换为等价方程x=ϕ(x),这样原方程的根必满足: x*=ϕ(x*),即ϕ(x)作用在x*上,其值不发生变化,因此我们也称x*为ϕ(x)的不动点,要求方程f(x)=0的根就转化为求ϕ(x)的不动点了。 具体迭代如下: 先取一个估计值x0来试探,若ϕ(x0)=x0,则x*=x0(可能性很小) 一般ϕ(x0)≠x0,记x1=ϕ(x0),若x1=ϕ(x1),则x*=x1 若x1≠ϕ(x1),记x2=ϕ(x1),再用x2继续试探…… 如此反复计算,即形成一迭代公式xk+1=ϕ(xk),(k=0,1,2,…) 如果{xk}收敛,则称迭代公式是收敛的;否则称迭代公式是发散的。 如果{xk}收敛于x*,而ϕ(x)是连续函数时,那么x*即是ϕ(x)的不动点。 也即x*就是方程的根。 图2 开始 输入 a,b,e y1= f(a),y2=f(b) y1*y2>0 输出失败信息 n=1 y1*y<0 b=cy2=y a=c,y1=y b-a n=n+1 输出x stop c=(a+b)/2,y=f(c) 是 是 否 否 是 否 5、Newton迭代法的思想: 给定方程f(x)=0,迭代公式为: ,应用该公式来解方程的方法就称为牛顿迭代法。 6、Newton迭代法的算法: 1)给出初始近似根x0及精度ε; 2)计算 →x1; 3)若|x1-x0|<ε或迭代次数达到预定指定的次数N,则转向4);否则x1→x0,转向2); 4)输出满足精度的根x1,结束。 (程序框图略) 7、收敛性判定定理: 定理1: 假定函数ϕ(x)满足下列条件: 1、对任意x∈[a,b]有ϕ(x)∈[a,b] 2、存在正数L<1,使对任意x1,x2∈[a,b]有 |ϕ(x1)-ϕ(x2)|≤L|x1-x2|(0≤L<1) 则迭代过程xk+1=ϕ(xk)对于任意初值x0∈[a,b]均收敛于方程x=ϕ(x)的根x*,且有如下的误差估计式: 局部收敛性定义: 如果存在不动点x*的某个邻域U(x*,δ),使得对于任意初值x0∈U(x*,δ),迭代公式xk+1=ϕ(xk)(k=0,1,2......)产生的数列{xk}均收敛与x*,则称迭代公式xk+1=ϕ(xk)是局部收敛的。 定理2: 设ϕ(x)在x=ϕ(x)的根x*邻近有连续的一阶导数,且|ϕ’(x*)|<1, 则迭代公式xk+1=ϕ(xk)具有局部收敛性。 8、收敛速度定义 定义当k→∞时,有 (C≠0,且为常数) 则称迭代过程是p阶收敛的。 特别地,当p=1,0 p=2称为平方收敛 其中: ek=(x*-xk)为迭代误差。 定理3: 对于迭代过程xk+1=ϕ(xk),如果ϕ(p)(x)在所求根x*的邻近连续,并且ϕ’(x*)=ϕ’’(x*)=...=ϕ(p-1)(x*)=0,ϕ(p)(x*)≠0,则该迭代过程在点x*邻近是P阶收敛的。 注: 牛顿迭代公式在根x*的邻近至少是平方收敛的。 三、数值计算实验(以下实验都需利用Matlab软件来完成) 实验2.1(求非线性方程根实验 (一)) 实验目的: 会求非线性方程的根 实验内容: 1、用Matlab软件做出下列方程的图形,观察根所在的区间: 2、用二分法求上述方程的根,并分析其误差。 实验2.1(求非线性方程根实验 (二)) 实验目的: 会求非线性方程的根 实验内容: 试有方程: 1、分别用Newton迭代法和快速弦截法求上述方程的根; 2、比较两种迭代法的收敛速度。 四、具体操作见下例。 (以下实验都需利用Matlab软件来完成) 例1: 用二分法求方程x2-x-1=0的正根,要求误差小于0.05 解: 1)首先根据程序框图编制二分法的函数文件: erfen.m 2)再编写函数文件fc.m funtiony=fc(x) y=x^2-x-1; 3)最后在Matlab命令窗口输入调用命令: erfen(‘fc’,0,10,0.05) 4)输出结果为n=56 ans=1.6180 例2: 求解非线性方程5x2sinx-e-x=0,观察知有多解,如何求出所有解? 解: 1)编写M—文件: newton.m 2)编写函数文件: fc1.m和df.m %fc1.m functiony=fc1(x) y=5*x^2*sin(x)-exp(-x); %df.m functiony=df(x) y=10*x*sin(x)+5*x^2*cosx+exp(-x); 3)用命令fplot(‘[5*x^2*sin(x)-exp(-x),0]’,[0,10]))作图,%在指定的范围内画出函数的图形。 (注意: [5*x^2*sin(x)-exp(-x),0]中的‘[…,0]’是作y=0直线,即x轴)(见图3) 从图中可看出方程在[0,10]区间有4个解,分别在0、3、6、9附近,所以 4)用命令 >>>>newton(x0) 即可求出求方程fc1=0在x0附近的近似解 得出结果 ans= 0.50183.14076.28329.4248 图3 注: 用Newton迭代法可求出多个根。 具体做法: 先用fplot命令作出函数的图形,再根据图形给出不同的初始值x0,最后通过调用Newton迭代法程序求出非线性方程的所有根。 五、思考题 思考题1、在二分法中取区间中点,每次计算一次函数值。 如果每次把区间三等份,计算两个函数值,是否可以找出方程的根? 给出它的算法。 和二分法比较它的效率如何? 思考题2、目的: 找出一维搜索的最佳途径。 内容: 假设f(x)=0在[a,b]区间内只有一个根(可以是重根),求解该方程等价于在[a,b]区间找|f(x)|的极小值点。 设计一种寻找极小值点的方法,使得计算f(x)的次数尽可能的少,并完成数值实验。 你能从理论上证明你的搜索方法最好吗? 思考题3、目的: 了解牛顿收敛法的结构和“局部”收敛性。 内容: 牛顿法可以直接用来求解复数方程z3-1=0,在复平面上它的三个根是z1*=1, 。 选择中心位于坐标原点边长为2的正方形内的点为初始值,把收敛到三个不同根的初始值分别标上不同的颜色。 只要计算足够多的点,你将得到关于牛顿收敛域的彩色图片。 第三章线性方程组的解法 一、主要要求 编写方程组求解的通用程序: Jacobi迭代、Seidel迭代以及SOR迭代程序 二、主要结果回顾 1、迭代法: •它的基本思想是将线性方程组AX=b化为: X=BX+f •再由此构造向量序列{X(k)}: X(k+1)=BX(k)+f •若{X(k)}收敛至某个向量x*,则可得向量X*就是所求方程组AX=b的准确解. •线性方程组的迭代法主要有Jocobi迭代法、Seidel迭代法和超松弛(Sor)迭代法. 2、收敛性判定定理: 定理1: 对任意初始向量X(0)及任意向量f,由此产生的迭代向量序列{X(k)}收敛的充要条件是 。 定理2: 若||B||<1,则迭代公式X(k)=BX(k-1)+f收敛. 3、Jacobi迭代公式: 其中: BJ=-D-1(U+L),fJ=D-1b。 4、Seidel迭代公式: 其中: Bs=-(D+L)-1U,fs=(D+L)-1b。 5、Jacobi迭代、Seidel迭代的算法: Jacobi迭代算法: 1)给出矩阵A,b,x0 2)计算B=D-1(-U-L),f=D-1b 3)y=BX0+f 4)若||y-x0||<=1.0e-6,则转到5),否则,转到3) 5)输出y和n。 Seidel迭代算法: 1)给出矩阵A,b,x0 2)计算B=(D-L)-1U,f=(D-L)-1b,则有: 3)y=BX0+f 4)若||y-x0||<=1.0e-6,则转到5),否则,转到3) 5)输出y和n。 6、程序框图: (见图4) 开始 Y N N 输入a,b,x0 由矩阵a,生成D,U,L,B,f y=B*x0+f n=1 ||y-x0||>=1.0e-6 x0=yy=B*x0+f n=n+1 输出x 结束 图4 仿Jacobi迭代法算法和Seidel迭代算法可给出超松驰迭代的算法(略) 三、数值计算实验(以下实验都需用Matlab软件来完成) 实验3.1(求解线性方程组实验) 实验目的: 掌握各种迭代法,比较各种迭代法的渐进收敛速度. 实验内容: 令 1、计算A的条件数cond(A);(可选用任何一种范数) 2、上述方程组是否为病态方程组? 若是,如何改变方程组的病态性? 3、分别用Jacobi迭代、Seidel迭代与超松驰迭代求方程组AX=b的近似解; 4、比较Jacobi迭代、Seidel迭代与超松驰迭代的渐进收敛速度。 注: 上述实验分两次完成。 四、具体操作见下例(以下实验都需用Matlab软件来完成) 例: 用Jacobi方法求解下列方程组,设X(0)=0,精度为10-6 解: 1)先根据Jacobi算法编写M—文件: Jacobi(a,b,X0) 2)在Matlab命令窗口输入命令: >>a=[10-10;-110-2;0-110]; >>b=[9;7;6]; >>Jacobi(a,b,[0;0;0]) 3)输出结果: y= 0.9938 0.9381 0.6938 n=9 注: n=9为迭代次数。 五、思考题 思考题: 目的: 以Hilbert矩阵为例,研究处理病态问题可能遇到的困难。 内容: 设Hilbert矩阵为 它是一个对称正定矩阵,而且cond(Hn)随着n的增加迅速增加。 其逆矩阵Hn-1=(aij),这里 1)画出ln(cond(Hn))~n之间的曲线(可以用任何一种范数)。 你能猜出ln(cond(Hn))~n之间有何种关系吗? 提出你的猜测并想法验证。 2)设D是Hn的对角元素开方构成的矩阵。 ,不难看出它仍然是对称矩阵,而且对角线元素都是1。 把Hn变成 的技术称为预处理(preconditioning)。 画出ln(cond( )/cond(Hn))~n之间的曲线(可以用任何一种范数)。 对于预处理你能得出什么印象? 3)对于4≤n≤12,给出不同的右端项b。 分别用X1=Hn-1b, ,X2=D-1 以及X3=Hn\b(这是Matlab的一条命令)求解HnX=b,比较计算结果。 4)取不同的n并以Hn的第一列为右端向量b,同高斯—塞德尔迭代求解HnX=b,观察其收敛性。 最后你能对于有关Hilbert矩阵的计算得出哪些结论。 第四章插值与拟合 一、主要要求 编写拉格朗日插值法、分段线性插值法、*三次样条插值通用程序(Matlab自带)。 拉格朗日插值多项式: 二、主要结果回顾 插值法是函数逼近的重要方法之一,有着广泛的应用。 在生产和实验中,函数f(x)或者其表达式复杂不便于计算或者无表达式而只有函数在给定点的函数值(或其导数值),此时我们希望建立一个简单的而便于计算的函数ϕ(x),使其近似的代替f(x),这就是所谓的插值法.有很多种插值法,其中以拉格朗日(Lagrange)插值和牛顿(Newton)插值为代表的多项式插值最有特点,常用的插值还有Hermit插值,分段插值和样条插值. 1、插值: 求近似函数的方法: 由实验或测量的方法得到所求函数y=f(x)在互异点x0,x1,...,xn处的值y0,y1,…,yn, 构造一个简单函数ϕ(x)作为函数y=f(x)的近似表达式: y=f(x)≈ϕ(x) 使ϕ(x0)=y0,ϕ(x1)=y1,⋯,ϕ(xn)=yn,(*) 这类问题称为插值问题。 f(x)称为被插值函数,ϕ(x)称为插值函数,x0,x1,...,xn称为插值节点。 (*)式称为插值条件。 2、Lagrange插值公式 其中x0,x1,...,xn为插值节点,yj为插值节点xj处的函数值(j=1,2,…n)。 3、Lagrange插值的截断误差(插值余项) 定理: 设Ln(x)是过点x0,x1,x2,…xn的n次插值多项式,f(n)(x)在[a,b]上连续,f(n+1)(x)在[a,b]上存在,其中[a,b]是包含点x0,x1,x2,…,xn的任一区间,则对任意给定的x∈[a,b],总存在一点∈ξ(a,b)(依赖于x)使 其中: 。 4、拉格朗日插值法程序框图: (见图5) (图5) 三、数值计算实验(以下实验都需利用Matlab软件来完成) 实验4.1(观察龙格(Runge)现象实验) 实验目的: 观察拉格朗日插值的龙格(Runge)现象. 实验内容: 对于函数 进行拉格朗日插值,取不同的节点数,在区间[-5,5]上取等距间隔的节点为插值点,把f(x)和插值多项式的曲线画在同一张图上进行比较。 (a可以取任意值)具体步骤: 1、a=1时,1)n=4,作出f(x)和插值多项式的曲线图; 2)n=10,作出f(x)和插值多项式的曲线图; 2、a=0.5时,1)n=4,作出f(x)和插值多项式的曲线图; 2)n=10,作出f(x)和插值多项式的曲线图; 3、分析上述曲线图,你可以得出什么结论? 四、具体操作 例1、给出f(x)=lnx的数值表,用Lagrange插值计算ln(0.54)的近似值。 x 0.4 0.5 0.6 0.7 0.8 Lnx -0.916291 -0.693147 -0.510826 -0.357765 -0.223144 解: 1)首先根据程序框图拉格朗日插值法的函数文件: lagrange(x0,y0,x) 2)在Matlab命令窗口输入调用命令: >>x0=[0.4: 0.1: 0.8]; >>y0=[-0.916291-0.693147-0.510826-0.356675-0.223144]; >>lagrange(x0,y0,0.54) 3)输出结果为: -0.616143(精确解: -0.616186) 实验4.2分段插值实验 实验目的: 分段插值计算,会使用一维插值函数,寻找最佳的插值方法。 实验内容: 在1—12的11小时内,每隔1小时测量一次温度,测得的温度依次为: 5,8,9,15,25,29,31,30,22,25,27,24。 1)试估计每隔1/2小时的温度值。 2)分别用分段线性插值,立方插值,三次样条插值和最邻近插值估计其值。 实验4.3(插值与拟合实验) 实验目的: 求下列数据的多项式拟合函数 x=[0,1,2,3,4,5,6,7,8,9,10]; y=[1,2.3,2.1,2,4.6,4.7,4.3,8.1,9.2,9.8,10.3]; 实验内容: 1、做出原始数据的散点图; 2、做出原始数据的折线图(即分段线性插值函数图); 3、做出原始数据的分段二次拟合曲线图(见图6); 4、做出原始数据的分段三次拟合曲线图。 (图6) 五、思考题 思考题1、目的: 观察最小二乘多项式的数值不稳定性现象。 内容: 1、在[-1,1]区间上取n=20个等距节点,计算出以相应节点上ex的值做为数据样本, 2、以1,x,x2,...,xt为基函数做l=3,5,7,9次的最小二乘多项式。 3、画出ln(cond(A))~n之间的曲线,其中A是确定最小二乘多项式系数的矩阵。 4、计算出不同阶最小二乘多项式给出的最小偏差 思考题2、目的: 观察对于非光滑函数进行多项式插值的可能性。 内容: 在[0,1]上取f(x)=|sinkπx|,选择不同的k和n,用等距节点做拉格朗日多插值项式,观察误差大小和收敛情况。 第五章数值积分 一、主要要求 编写定步长复合梯形公式、定步长复合Simpson公式、变步长梯形法及龙贝格算法*通用子程序; 二、主要结果回顾 1、梯形公式: 对于连续函数f(x),有积分中值定理 其中f(ξ)是被积函数f(x)在积分区间上的平均值。 因此,如果我们能给出求平均值f(ξ)的一种近似方法,相应地就可以得到计算定积分的一种数值方法。 如果取 ,即得下列梯形公式: 二、辛普生(Simpson)公式: 先用某个简单函数近似逼近f(x),然后用 在[a,b]区间的积分值近似表示f(x)在[a,b]区间上的定积分,即取 我们可以取 为前面介绍的f(x)的插值多项式Ln(x)或拟合多项式P(x)进行近似计算。 如取 为插值多项式Ln(x),则相应得到的数值积分公式称为插值型求积公式。 如: 若考虑过点A,B,C三点的抛物线段L代替曲线段y=f(x)(见图7)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 实验 指导书