第一章 误差与范数

1.1 误差的来源例1.1.1 用差商ha f h a f a f )()()(-+≈'求x x f ln )(=在3=x 处导数的近似值.取1.0=h ，1000.0=h ，h =0.000 000 000 000 001和h =0.000 000 000 000 000 1分别用MATLAB 软件计算，取十五位数字计算.解 在MATLAB 工作窗口输入下面程序>>a=3;h=0.1;y=log(a+h)-log(a);yx=y/h运行后得yx = 0.32789822822991 将此程序中h 改为0.000 1，运行后得yx = 0.33332777790385后者比前者好.再取h = 0.000 000 000 000 001，运行后得yx = 0.44408920985006不如前者好.取h = 0.000 000 000 000 000 1，运行后得yx = 0算出的结果反而毫无价值.例1.1.2 分别求方程组b AX =在下列情况时的解，其中A ⎪⎪⎭⎫⎝⎛=011111.. （1）⎪⎪⎭⎫ ⎝⎛=22b ;（2）⎪⎪⎭⎫⎝⎛=0122.b . 解 （1） 首先将方程组b AX =化为同解方程b A X 1-=，然后在MATLAB 工作窗口输入程序>> b=[2,2]';A=[1,1;1,1.01]; X=A\b运行后输出当⎪⎪⎭⎫ ⎝⎛=22b 时，b AX =的解为⎪⎪⎭⎫⎝⎛=02X ；（2）同理可得，当⎪⎪⎭⎫ ⎝⎛=0122.b 时，b AX =的解为⎪⎪⎭⎫⎝⎛=11X .例1.1.3 计算e 的近似值. 解 泰勒级数e +++++++=!!4!3!21432 n x x x x x n x)(∞<<-∞x , 取1=x ，得e +++++++=!1!41!31!2111　　n . （1.2）这是一个无限过程，计算机无法求到精确值.只能在（1.2）取有限项时计算，再估计误差.如果取有限项!!!!)(　　n s n 1413121111++++++=作为e 的值必然会有误差，根据泰勒余项定理可知其截断误差为e !)1()1(　+=-n e s n θ)10(<<θ.如果取（1.2）的前九项，输入程序>> n =8;s=1;S =1; for k=1:ns=s*k;S=S+1/s, ends,S,R=3/(s*(n+1)) 或>>S1=1+1+1/2+1/(1*2*3)+1/(1*2*3*4)+1/(1*2*3*4*5)+1/(1*2*3*4*5*6)+1/(1*2*3*4*5*6*7)+1/(1*2*3*4*5*6*7*8),R1=3/(1*2*3*4*5*6*7*8*9)运行后结果S = R =2.71827876984127 8.267195767195768e-006 因为截断误差为e ),10(101968.267!93!)18()1(6-8<<⨯≈<+=-θθ　　e s 所以e 的近似值e ≈≈++++++++=!81!71!61!51!41!31!2111)1(8 s 2.718 28.1.2 误差和有效数字例1.2.1 取282.718作为e 的四舍五入近似值时，求其绝对误差和相对误差. 解 在MATLAB 工作窗口输入程序>>juewu=exp(1)-2.71828运行后输出结果为juewu = 1.828 459 045 505 326e-006例1.2.2 计算⎰π20sin x xd x 的近似值，并确定其绝对误差和相对误差.解 因为被积函数xxsin 的原函数不是初等函数，故用泰勒级数求之.++-+-=!!!!sin ９　７５　３８６x x x x x x 421 )(∞<<-∞x , （1.5） 这是一个无限过程，计算机无法求到精确值.可用（1.5）的前四项!!!７５　３６x x x -+-421代替被积函数xxsin ，得 ⎰π=20sin x x y d ⎰π≈20(x !!!142７５　３６x x x -+-)d x =!7)2(!5)2(!3)2(2753７５　３⋅π-⋅π+⋅π-π=y ˆ. 根据泰勒余项定理和交错级数收敛性的判别定理，得到绝对误差!99)2(ˆ9⋅<-=πyy R = WU ， 在MA TLAB 命令窗口输入计算程序如下：>>syms xf=1-x^2/(1*2*3)+x^4/(1*2*3*4*5)-x^6/(1*2*3*4*5*6*7)y=int(f,x,0,pi/2),y1=double(y)y11=pi/2-(pi/2)^3/(3*3*2)+(pi/2)^5/(5*5*4*3*2)-(pi/2)^7/(7*7*6*5*4*3*2)inf=int(sin(x)/x,x,0,pi/2) ,infd=double(inf) WU =(pi/2)^9/(9*9*8*7*6*5*4*3*2), R =infd-y11因为运行后输出结果为： =y 1.370 762 168 154 49，yˆ=1.370 744 664 189 38，=R 1.750 396 510 491 47e-005， WU = 1.782 679 830 970 664e-005410-<.所以，yˆ的绝对误差为=ε410-，故⎰π=20sin x xy d 7 1.370≈x .yˆ的相对误差为 =r ε71.37010ˆ4-=y ε<0.007 3%.1.3 误差估计的基本方法例1.3.4 设计三种算法求方程01522=-+x x 在)3,2(的一个正根*x 的近似值，并研究每种算法的误差传播情况.解 为解已知方程，我们可以设计如下三种算法，然后将计算结果与此方程的精确解5.2*=x 比较，观察误差的传播.算法1 将已知方程化为同解方程=x 2215x -.取初值20=x ，按迭代公式21215k k x x -=+依次计算 ,,,,21n x x x ，结果列入表1–3中.算法2 将已知方程化为同解方程1215+=x x .取初值20=x ，按迭代公式 12151+=+k k x x依次计算 ,,,,21n x x x ，结果列入表1–3中.算法3 将已知方程化为同解方程141522+-+-=x x x x x .取初值20=x ，按迭代公式为1415221+-+-=+k k kk k x x x x x 依次计算 ,,,,21n x x x ，结果列入表1–3中.我们为这三种算法的计算编写两套MATLAB 程序如下： （1）MATLAB 主程序function [k,juecha,xiangcha,xk]= liti112(x0,x1,limax) % 输入的量--x0是初值, limax 是迭代次数和精确值x; % 输出的量--每次迭代次数k 和迭代值xk,% --每次迭代的绝对误差juecha 和相对误差xiangcha ， x(1)=x0;for i=1:limaxx(i+1)=fl(x(i));%程序中调用的fl.m juecha = abs(x(i)-x1);xiangcha = juecha /( abs(x(i))+eps);xk=x(i);k=i-1;[(i-1),juecha,xiangcha,xk] end（2）MATLAB 调用函数程序及其计算结果①算法2的MATLAB 调用函数程序function y1=fl(x)y1=15/(2*x+1);② 将MATLAB 主程序和调用函数程序分别命名liti112.m 和fl.m ，分别保存为M 文件，然后在MATLAB 工作窗口输入命令>> [k,juecha,xiangcha,xk]= liti112(2,2.5,100) ③运行后输出计算结果列入表1–3和表 1-4中.④将算法2的MATLAB 调用函数程序的函数分别用y1=15-2*x^2和y1=x-(2*x^2+x-15)/(4*x+1)代替，得到算法1和算法3的调用函数程序，将其保存，运行后将三种算法的前8个迭代值821,,,x x x 列在一起（见表 1-3），进行比较.将三种算法的821,,,x x x 对应的绝对误差和相对误差的值列在一起（见表 1-4），进行比较.1.4 数值计算中应注意的问题例1.4.1 求数）181(71915-+⨯=-x 的近似值. 解 （1）直接用MATLAB 命令>> x=(7^15)*(sqrt(1+8^(-19))-1)运行后输出结果x = 0问题出现在两个相近的数1981-+与1相减时，计算机运行程序>>sqrt(1+8^(-19))-1运行后输出结果ans = 0由于计算机硬件只支持有限位机器数的运算，因此在计算中可能引入和传播舍入误差.因为有效数字的严重损失，导致输出18119-+-的结果为0，计算机不能再与数157继续进行真实的计算，所以，最后输出的结果与x 的精确值不符.（2）如果化为18187)181(71919151915++⨯=-+⨯=---x ，再用MATLAB 命令>> x=(7^15)*( (8^(-19))/(sqrt(1+8^(-19))+1))运行后输出结果x = 1.6471e-005 这是因为18119-+-化为18181919++--后，计算机运行程序>> x= (8^(-19))/(sqrt(1+8^(-19))+1)运行后的结果为x =3.4694e-018 由于有效数字的损失甚少，所以运算的结果-18103.4694⨯再与157继续计算，最后输出的结果与x 的精确值相差无几.例1.4.2 求数)13030ln(2--=y 的近似值. 解 （1）直接用MATLAB 程序>> x=30;x1= sqrt(x^2-1)运行后输出结果x1 = 29.9833 输入MATLAB 程序>> x=30; x1=29.9833;y=log(x-x1)运行后输出结果y = -4.0923（2）因为)13030ln(2--中的30=x 很大，如果采用倒数变换法111221-+=--=x x x x z ，即130301ln)13030ln(22-+=--)190030ln(-+-=.输入MATLAB 程序>> x=30;y=-log(x+sqrt(x^2-1))运行后输出结果y = -4.0941（3）输入MA TLAB 程序>> x=30; y=log(x-sqrt(x^2-1))运行后输出结果y = -4.0941 可见，（2）计算的近似值比（1）的误差小.参加计算的数，有时数量级相差很大.如果不注意采取相应的措施，在它们的加减法运算中，绝对值很小的那个数经常会被绝对值较大的那个数“吃掉”，不能发挥其作用，造成计算结果失真.例1.4.4 请在16位十进制数值精度计算机上利用软件MATLAB 计算下面的两个数0.30.1111111111111111*++=x 和0.30.11111111111111111*++=y将计算结果与准确值比较，解释计算结果.解 在MATLAB 工作窗口输入下面程序>> x=111111*********+0.1+0.3, y=1111111111111111+0.1+0.3运行后输出结果x = 1.111111*********e+014，y =1.111111*********e+015 从输出的结果可以看出，x *x =，而y *y ≠.为什么*y 仅仅比*x 多一位1，而y *y ≠呢？这是因为计算机进行运算时，首先要把参加运算的数写成绝对值小于1而“阶码”相同的数，这一过程称为数的“对阶”.在16位十进制数值精度计算机上利用软件MATLAB 计算这两个数，把运算的数*x 写成浮点规格化形式为，151515*103000**********.0001010000000000000.000100111111111111111.0⨯+⨯+⨯=x在16位十进制数值精度计算机上，三项的数都表示为小数点后面16位数字的数与1510之积，所以，计算机没有对数进行截断，而是按原来的三个数进行计算.因此，计算的结果x *x =.而161616*10030000000000000.00010010000000000000.000101111111111111111.0⨯+⨯+⨯=y三项的数都表示写成绝对值小于1而“阶码”都为1610的数以后，第一项的纯小数的小数点后面有16位数字.但是，后两项数的纯小数的小数点后面有17位数字，超过了16位十进制数值精度计算机的存储量，计算机对后两项的数都进行截断最后一位，即后两项的数都是16位机上的零，再进行计算，所以计算结果与实际不符.五、向量和矩阵的范数例1．5.1: 用matlab 求下列向量的2,1,5，,-∞∞范数。

习 题1.11. 分别求方程组b AX=在下列情况时的解，其中⎪⎪⎭⎫⎝⎛=01.2222A .（1）⎪⎪⎭⎫ ⎝⎛=22b ；（2）⎪⎪⎭⎫ ⎝⎛=01.22b . 2. 用差商h a f h a f a f )()()('-+≈，求x x f ln )(=在5=x 处导数)5('f 的近似值，取十五位数字计算，h 分别取0.000 1和0.000 000 01.习 题1.21. 为了使 4.420=的近似值的相对误差小于%1.0，试问取几位有效数字？2. 写出光速10*10)001000.0925997.2(⨯±=c s /cm 的绝对误差和相对误差. 3. 取655923.141作为π的四舍五入近似值时，求其有效数字、绝对误差和相对误差.4. 问722,141.3,142.3分别作为π的近似值各具有几位有效数字？习 题1.31.测得电压)2110(±=V V ，电流)5.020(±=I A ，则由欧姆定律得Ω==5.5I VR ，求R 的绝对误差和相对误差.2. 设nx y =，求y 的相对误差与x 的相对误差的关系式.3. 设计三种算法求方程0322=-+x x 在)2,1(的一个正根*x 的近似值，并研究每种算法的误差传播情况.4. 已测得某场地长*x 的值110=x m ，宽*y 的值80=y m ，已知x 的绝对误差是2.0m ，y 的绝对误差是1.0m ，求场地面积xy S =的绝对误差和相对误差.5. 正方形的边长约为200cm ，怎样测量才能使其面积误差不超过2cm 1？1. 求数）274(41915-+⨯=-x 的近似值.2. 利用四位数学用表求2cos 1-的近似值，要求至少有2位有效数字.3. 请在16位十进制数值精度计算机上利用软件MATLAB 计算下面的两个数0.30.1333333333333333*++=x 和0.30.133333333333333333*++=y将计算结果与准确值比较，解释计算结果.4. 利用恒等变换，使下列表达式的计算结果比较精确.（1）1,11>>--+x x x x x ； （2）e12-x ； （3）x xsin cos 1-. 5. 设x y ln =.当)0(>≈a a x 时，如果已知对数a ln 的绝对误差限为n-⨯1021，试估计真数a 的相对误差限及有效数字位数.6. 指出下列各题的合理计算途径（对给出具体数据的，请算出结果）：（1）o1cos 1-（三角函数值取四位有效数字）；（2）)13030ln(2--；（3）x xsin cos 1- （其中x充分小）； （4）127x .习 题 1.51.设⎪⎪⎪⎪⎪⎭⎫⎝⎛=1095791068565778710A ,当b T)31,33,23,32(=有微小误差)1.0,1.0,1.0,1.0(--=δb T 时，估计方程组b Ax =解的变化.2. 用MATLAB 求矩阵A 的2范数条件数、1范数条件数、∞范数条件数、Frobenius 范数条件数、条件数倒数的估计值,其中3. 用MATLAB 求第2题中矩阵A 的2范数、1范数、∞范数、Frobenius 范数、谱半径.4. 用MATLAB 求下列向量的2范数、1范数、∞范数、5范数、-∞范数:(1))5,2,1,1(-=X ;(2) )12,45,06.0,89,4.0(--=Y T. 5. 线性方程组求解与性态讨论.求b Ax =的解向量x ，其中然后把b 扰动为b ˆT )99.30,01.33,99.22,01.32(=，再求解b xA ˆˆ=. 计算xx x -ˆ（使用1范数或∞范数），讨论方程组性态..98755.03.02.064321543210⎪⎪⎪⎪⎪⎭⎫⎝⎛---=A .31332332,9105791068565778710⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫⎝⎛=b A。

第一章误差与范数1.1 误差的来源1.2 误差和有效数字1.3 误差估计的基本方法1.4 数值计算中应注意的问题第二章非线性方程（组）的数值解法2.1 方程（组）的根及其MATLAB命令2.1.2 求解方程（组）的solve命令2.1.4 求解方程（组）的fsolve命令2.2 搜索根的方法及其MATLAB程序2.2.1 作图法及其MATLAB程序2.2.2 逐步搜索法及其MATLAB程序2.3 二分法及其MATLAB程序2.3.1 二分法求方程根的迭代次数的MATLAB命令2.3.2 二分法的MATLAB程序2.4 迭代法及其MATLAB程序2.4.2 迭代法的MATLAB程序12.4.5 迭代法的MATLAB程序22.5 迭代过程的加速方法及其MATLAB程序2.5.2 加权迭代的MATLAB程序2.5.4 艾特肯(Aitken)加速方法的MATLAB程序2.6 牛顿(Newton)切线法及其MATLAB程序2.6.2 牛顿切线法的收敛性及其MATLAB程序2.6.3 牛顿切线法的MATLAB程序2.6.6 牛顿切线法的加速及其两种MATLAB程序2.7 割线法及其MATLAB程序2.7.2 割线法的MATLAB程序2.8 抛物线法及其MATLAB程序2.8.2 抛物线法的MATLAB程序2.9 求解非线性方程组的牛顿法及其MATLAB程序2.9.1 求解二阶非线性方程组的牛顿法及其MATLAB程序2.9.2 求解n阶非线性方程组的牛顿法及其MATLAB程序2.9.3 求解n阶非线性方程组的拟牛顿法及其MATLAB程序第三章解线性方程组的直接方法3.1 方程组的逆矩阵解法及其MATLAB程序3.1.3 线性方程组有解的判定条件及其MATLAB程序3.2 三角形方程组的解法及其MATLAB程序3.2.2 解三角形方程组的MATLAB程序3.3 高斯消元法和列主元素消元法及其MATLAB程序3.3.1 高斯消元法及其MATLAB程序3.3.2 列主元素消元法及其MATLAB程序3.4 LU分解法及其MATLAB程序3.4.1 判断矩阵LU分解的充要条件及其MATLAB程序3.4.2 直接LU分解法及其MATLAB程序3.4.3 含交换阵的LU分解法及其MATLAB程序3.4.4 正定对称矩阵的Cholesky分解及其MATLAB程序3.4.5 正定对称矩阵的Cholesky分解及其MATLAB程序3.5 求解线性方程组的LU方法及其MATLAB实现3.5.1 求解线性方程组的Cholesky分解及其MATLAB实现3.5.2 求解线性方程组的直接LU分解法及其MATLAB实现3.5.3 解线性方程组的选主元的LU方法及其MATLAB实现3.6 误差分析及其两种MATLAB程序3.6.1 用MATLAB软件作误差分析AX=解的性态的MATLAB程序3.6.2 求P条件数和讨论bAX=解和A的性态的MATLAB程序3.6.3 用P范数讨论b第四章解线性方程组的迭代法4.1 迭代法和敛散性及其MATLAB程序4.1.2 用谱半径判别迭代法产生的迭代序列的敛散性的MATLAB程序4.2 雅可比（Jacobi）迭代及其MATLAB程序4.2.2 判别雅可比迭代收敛性的MATLAB程序4.2.3 雅可比迭代的两种MATLAB程序4.3 高斯—塞德尔（Gauss-Seidel）迭代及其MATLAB程序4.3.3 高斯—塞德尔迭代两种MATLAB程序4.4 解方程组的超松弛迭代法及其MATLAB程序4.4.2 超松弛迭代法收敛性及其MATLAB程序4.4.3 超松弛迭代法的MATLAB程序第五章矩阵的特征值与特征向量的计算5.2 幂法及其MATLAB程序5.2.2 幂法的MATLAB程序5.3 反幂法和位移反幂法及其MATLAB程序5.3.3 原点位移反幂法的两种MATLAB程序5.4 雅可比(Jacobi)方法及其MATLAB程序5.4.3 雅可比方法的MATLAB程序5.5 豪斯荷尔德方法及其MATLAB程序5.5.1 豪斯荷尔德方法及其MATLAB程序5.5.2 矩阵约化为上Householder阵及其MATLAB程序5.5.3 实对称矩阵的三对角化及其MATLAB 程序5.6 QR 方法及其MATLAB 程序5.6.2 QR 方法的MATLANB 程序5.6.5 最末元素位移QR 法计算实对称矩阵特征值及其MATLAB 程序5.6.6 求根位移QR 方法计算实对称矩阵A 的特征值及其MATLAB 程序5.7 广义特征值问题及其MATLAB 程序5.7.2 用MATLAB 计算BX AX λ=型的广义特征值和特征向量第六章 函数的插值方法6.1 插值问题及其误差6.1.2 与插值有关的MATLAB 函数6.2 拉格朗日（Lagrange ）插值及其MATLAB 程序6.2.1 线性插值及其MATLAB 程序6.2.2 抛物线插值及其MATLAB 程序 6.2.3 n 次拉格朗日（Lagrange ）插值及其MATLAB 程序6.2.5 拉格朗日多项式和基函数的MATLAB 程序6.2.6 拉格朗日插值及其误差估计的MATLAB 程序6.3 牛顿（Newton ）插值及其MATLAB 程序6.3.3 牛顿插值多项式、差商和误差公式的MATLAB 程序6.3.4 牛顿插值及其误差估计的MATLAB 程序6.3.5 牛顿插值法的MATLAB 综合程序6.4 埃尔米特（Hermite ）插值及其MATLAB 程序6.4.3 埃尔米特插值多项式的MATLAB 程序6.5 分段插值及其MATLAB 程序6.5.1 高次插值的振荡和MATLAB 程序6.5.3 分段线性插值的MATLAB 程序6.5.4 作有关分段线性插值图形的MATLAB 程序6.5.5 用MATLAB 计算有关分段线性插值的误差6.6 分段埃尔米特（Hermite ）插值及其MATLAB 程序6.6.2 Hermite 插值的MATLAB 程序6.6.3 作有关分段Hermite 插值图形的MATLAB 程序6.6.4 用MATLAB 计算有关分段Hermite 插值的误差6.7 三次样条插值及其MATLAB 程序6.7.4 用一阶导数计算的几种样条函数和MATLAB 程序6.7.6 用MATLAB 计算三次样条6.7.7 几种作三次样条有关图像的MATLAB 程序6.7.8 用MATLAB 计算有关分段三次样条的误差6.8 高元插值及其MATLAB 程序6.8.1 MESHGRID 命令的功能和调用格式6.8.2 单调数据点上的二元插值及其MATLAB 程序6.8.3 三元插值及其MATLAB 程序第七章 函数逼近与曲线（面）拟合7.1 曲线拟合、误差及其MATLAB 实现7.2 曲线拟合的线性最小二乘法及其MATLAB 程序7.3 函数)(x r k 的选取及其MATLAB 程序7.4 多项式拟合及其MATLAB 程序7.5 拟合曲线的线性变换及其MATLAB 程序7.6 函数逼近及其MATLAB 程序7.7 三角多项式逼近及其MATLAB 程序7.8 随机数据点上的二元拟合及其MATLAB 程序7.9 随机数据点上的n 元拟合及其MATLAB 程第八章 数值微分8.2 一阶导数的数值计算及其MATLAB 程序 (117)8.2.1 差商求导及其MATLAB 程序8.2.2 中心差商公式求导四种MATLAB 程序8.2.3 Richardson 外推法求导及其MATLAB 程序8.2.4 牛顿多项式求导及其MATLAB 程序8.2.5 diff 和gradient 函数在数值求导中的应用8.3 高阶导数的数值计算及其MATLAB 程序8.3.1 插值或拟合高阶数值导数及其MATLAB 程序8.3.2 高阶泰勒数值导数及其MATLAB 程序8.4 数值梯度和数值偏导数的计算及其MATLAB 程序8.4.1 梯度和偏导数的数值计算及其MATLAB 程序8.4.2 计算雅克比矩阵及其行列式的MATLAB 方法第九章 数值积分9.1 积分的符号计算9.1.1 定积分的MATLAB 符号计算9.1.2 变限积分的MATLAB 符号计算9.2 数值积分的思想及其MATLAB 程序9.2.3 矩形公式的MATLAB 运算9.3 插值型数值积分及其MATLAB 程序9.3.2 梯形公式的四种MATLAB 程序9.3.4 辛普生（Simpson ）数值积分的MATLAB 程序9.3.6 牛顿-科茨（Newton-Cotes ）数值积分和误差分析的三种MATLAB 程序9.3.7 利用三次样条求表格型数值积分的MATLAB 方法9.3.8 利用拉格朗日插值等方法求表格型数值积分的MATLAB 方法9.4 龙贝格公式及其MATLAB 程序9.4.2 龙贝格积分的MATLAB 程序9.5 自适应积分及其MATLAB 程序9.6 高斯（Gauss ）型积分公式及其MATLAB 程序9.6.2 在]1,1[ 上的高斯-- 勒让德积分公式及其MATLAB 程序9.6.3 在],[b a 上的高斯-- 勒让德积分公式及其MATLAB 程序9.6.4 Radau 和Lobatto 积分公式及其MATLAB 程序9.7 反常积分的计算及其MATLAB 程序9.7.1 无穷积分的符号计算及其MATLAB 程序9.7.2 无穷积分的近似计算及其MATLAB程序9.7.3 无界函数反常积分的符号计算的五种MATLAB程序9.7.4 无界函数反常积分的近似计算的两种MATLAB程序9.8 多重积分的计算及其MATLAB程序9.8.1 二重积分的符号计算及其MATLAB程序9.8.2 二重积分的梯形公式及其MATLAB程序9.8.3 矩形域上的辛普生公式及其MATLAB程序9.8.4 一般域上二重积分的数值计算及其MATLAB程序9.8.5 三重积分的计算及其MATLAB程序第十章常微分方程（组）求解10.1 常微分方程（组）的MATLAB符号求解10.1.1 MATLAB求常微分方程（组）的通解10.1.2 MATLAB求常微分方程（组）的特解10.1.3 线性常微分方程组的MATLAB解法10.3 欧拉方法及其MATLAB程序10.3.2 向前欧拉方法的三种MATLAB程序10.3.4 向后欧拉方法的MATLAB程序10.4 改进的欧拉方法及其MATLAB程序10.4.2 梯形公式的两个MATLAB程序10.4.4 改进的欧拉公式的MATLAB程序10.5 龙格—库塔方法10.5.2 二阶龙格—库塔方法及其MATLAB程序10.5.3 三阶龙格—库塔方法及其MATLAB程序10.5.4 四阶龙格—库塔方法及其MATLAB程序10.5.5 自适应龙格—库塔方法及其MATLAB库函数10.6 线性多步法及其MATLAB程序10.6.2 亚当斯显式公式及其MATLAB程序10.6.3 亚当斯隐式公式及其MATLAB程序10.6.4 米尔恩公式及其MATLAB程序10.6.5 汉明公式及其MATLAB程序10.6.6 预测-校正系统及其MATLAB程序10.7 一阶（高）阶微分方程（组）的数值解及其MATLAB程序10.7.1 一阶微分方程组的数值解及其MATLAB程序10.7.2 高阶微分方程（组）的数值解及其MATLAB实现10.8 边值问题的数值解及其MATLAB程序10.8.1 打靶法及其MATLAB程序10.8.2 有限差分方法及其MATLAB程序10.8.3 求解常微分方程（组）边值问题数值解的MATLAB库函数。

第一章 误差分析与向量与矩阵的范数一、内容提要本章要求掌握绝对误差、相对误差、有效数字、误差限的定义及其相互关系；掌握数值稳定性的概念、设计函数计算时的一些基本原则和误差分析；熟练掌握向量和矩阵范数的定义及其性质。

1．误差的基本概念和有效数字 1）．绝对误差和相对误差的基本概念设实数x 为某个精确值，a 为它的一个近似值，则称a x -为近似值a 的绝对误差，简称为误差． 当0≠x 时，xax -称为a 的相对误差．在实际运算中，精确值x 往往是未知的，所以常把aax -作为a 的相对误差． 2）．绝对误差界和相对误差界的基本概念设实数x 为某个精确值，a 为它的一个近似值，如果有常数a e ，使得 a e a x ≤- 称a e 为a 的绝对误差界，或简称为误差界．称ae a 是a 的相对误差界．此例计算中不难发现，绝对误差界和相对误差界并不是唯一的，但是它们越小，说明a 近似x 的程度越好，即a 的精度越好．3）．有效数字设实数x 为某个精确值，a 为它的一个近似值，写成n k a a a a 21.010⨯±=它可以是有限或无限小数的形式，其中),2,1( =i a i 是9,,1,0 中的一个数字，k a ,01≠为整数．如果n k a x -⨯≤-1021则称a 为x 的具有n 位有效数字的近似值．如果a 有n 位有效数字，则a 的相对误差界满足：n a a a x -⨯≤-111021。

4）．函数计算的误差估计如果),,,(21n x x x f y =为n 元函数，自变量n x x x ,,,21 的近似值分别为n a a a ,,,21 ，则)(),,,(),,,(12121k k n k akn n a x x fa a a f x x x f -⎪⎪⎭⎫⎝⎛∂∂≈-∑= 其中),,,(21n kak a a a f x x f ∂∂=⎪⎪⎭⎫ ⎝⎛∂∂，所以可以估计到函数值的误差界，近似地有 k a n k aka n n e x fe a a af x x x f ∑=⎪⎪⎭⎫⎝⎛∂∂≈≤-12121),,,(),,,( 如果令2=n ，设21,x x 的近似值分别为21,a a ，其误差界为111a e a x ≤-和≤-22a x 2a e ，取),(21x x f y =为21,x x 之间的四则运算，则它们的误差估计为，1121a a a a e e e +≈±；112121a a a a e a e a e +≈⋅；22211121a e a e a e a a a a +≈，02≠a 。

误差的来源例1.1.1 用差商ha f h a f a f )()()(-+≈'求x x f ln )(=在3=x 处导数的近似值.取1.0=h ，1000.0=h ，h =0.000 000 000 000 001和h =0.000 000 000 000 000 1分别用MATLAB 软件计算，取十五位数字计算.解 在MATLAB 工作窗口输入下面程序>>a=3;h=0.1;y=log(a+h)-log(a);yx=y/h运行后得yx = 0.32789822822991 将此程序中h 改为0.000 1，运行后得yx = 0.33332777790385后者比前者好.再取h = 0.000 000 000 000 001，运行后得yx = 0.44408920985006不如前者好.取h = 0.000 000 000 000 000 1，运行后得yx = 0算出的结果反而毫无价值.例1.1.2 分别求方程组b AX =在下列情况时的解，其中A ⎪⎪⎭⎫⎝⎛=011111.. （1）⎪⎪⎭⎫ ⎝⎛=22b ;（2）⎪⎪⎭⎫⎝⎛=0122.b . 解 （1） 首先将方程组b AX =化为同解方程b A X 1-=，然后在MATLAB 工作窗口输入程序>> b=[2,2]';A=[1,1;1,1.01]; X=A\b运行后输出当⎪⎪⎭⎫ ⎝⎛=22b 时，b AX =的解为⎪⎪⎭⎫⎝⎛=02X ；（2）同理可得，当⎪⎪⎭⎫ ⎝⎛=0122.b 时，b AX =的解为⎪⎪⎭⎫⎝⎛=11X .例1.1.3 计算e 的近似值. 解 泰勒级数e +++++++=!!4!3!21432 n x x x x x n x)(∞<<-∞x , 取1=x ，得e +++++++=!1!41!31!2111　　n . （1.2）这是一个无限过程，计算机无法求到精确值.只能在（1.2）取有限项时计算，再估计误差.如果取有限项!!!!)(　　n s n 1413121111++++++=作为e 的值必然会有误差，根据泰勒余项定理可知其截断误差为e !)1()1(　+=-n e s n θ)10(<<θ.如果取（1.2）的前九项，输入程序>> n=8; s=1;S =1; fork=1:n s=s*k; S=S+1/s, end s, S,R=3/(s*(n+1)) 或>>S1=1+1+1/2+1/(1*2*3)+1/(1*2*3*4)+1/(1*2*3*4*5)+1/(1*2*3*4*5*6)+1/(1*2*3*4*5*6*7)+1/(1*2*3*4*5*6*7*8),R1=3/(1*2*3*4*5*6*7*8*9)运行后结果S = R =2.71827876984127 8.267195767195768e-006 因为截断误差为e ),10(101968.267!93!)18()1(6-8<<⨯≈<+=-θθ　　e s 所以e 的近似值e ≈≈++++++++=!81!71!61!51!41!31!2111)1(8 s 2.718 28.1.2 误差和有效数字例1.2.1 取282.718作为e 的四舍五入近似值时，求其绝对误差和相对误差. 解 在MATLAB 工作窗口输入程序>>juewu=exp(1)-2.71828运行后输出结果为juewu = 1.828 459 045 505 326e-006例1.2.2 计算⎰π20sin x xd x 的近似值，并确定其绝对误差和相对误差.解 因为被积函数xxsin 的原函数不是初等函数，故用泰勒级数求之.++-+-=!!!!sin ９　７５　３８６x x x x x x 421 )(∞<<-∞x , （1.5） 这是一个无限过程，计算机无法求到精确值.可用（1.5）的前四项!!!７５　３６x x x -+-421代替被积函数xxsin ，得 ⎰π=20sin x x y d ⎰π≈20(x !!!142７５　３６x x x -+-)d x =!7)2(!5)2(!3)2(2753７５　３⋅π-⋅π+⋅π-π=y ˆ. 根据泰勒余项定理和交错级数收敛性的判别定理，得到绝对误差!99)2(ˆ9⋅<-=πyy R = WU ， 在MA TLAB 命令窗口输入计算程序如下：syms xf=1-x^2/(1*2*3)+x^4/(1*2*3*4*5)-x^6/(1*2*3*4*5*6*7) y=int(f,x,0,pi/2),y1=double(y)y11=pi/2-(pi/2)^3/(3*3*2)+(pi/2)^5/(5*5*4*3*2)-(pi/2)^7/(7*7*6*5*4*3*2)inf=int(sin(x)/x,x,0,pi/2) ,infd=double(inf) WU =(pi/2)^9/(9*9*8*7*6*5*4*3*2), R =infd-y11因为运行后输出结果为： =y 1.370 762 168 154 49，yˆ=1.370 744 664 189 38，=R 1.750 396 510 491 47e-005， WU = 1.782 679 830 970 664e-005410-<.所以，yˆ的绝对误差为=ε410-，故⎰π=20sin x xy d 7 1.370≈x .yˆ的相对误差为 =r ε71.37010ˆ4-=y ε<0.007 3%.1.3 误差估计的基本方法例1.3.4 设计三种算法求方程01522=-+x x 在)3,2(的一个正根*x 的近似值，并研究每种算法的误差传播情况.解 为解已知方程，我们可以设计如下三种算法，然后将计算结果与此方程的精确解5.2*=x 比较，观察误差的传播.算法1 将已知方程化为同解方程=x 2215x -.取初值20=x ，按迭代公式21215k k x x -=+依次计算 ,,,,21n x x x ，结果列入表1–3中.算法2 将已知方程化为同解方程1215+=x x .取初值20=x ，按迭代公式 12151+=+k k x x依次计算 ,,,,21n x x x ，结果列入表1–3中.算法3 将已知方程化为同解方程141522+-+-=x x x x x .取初值20=x ，按迭代公式为1415221+-+-=+k k kk k x x x x x 依次计算 ,,,,21n x x x ，结果列入表1–3中.我们为这三种算法的计算编写两套MATLAB 程序如下： （1）MATLAB 主程序function [k,juecha,xiangcha,xk]= liti112(x0,x1,limax) % 输入的量--x0是初值, limax 是迭代次数和精确值x; % 输出的量--每次迭代次数k 和迭代值xk,% --每次迭代的绝对误差juecha 和相对误差xiangcha ， x(1)=x0;for i=1:limaxx(i+1)=fl(x(i));%程序中调用的fl.m juecha = abs(x(i)-x1);xiangcha = juecha /( abs(x(i))+eps);xk=x(i);k=i-1;[(i-1),juecha,xiangcha,xk] end（2）MATLAB 调用函数程序及其计算结果 ①算法2的MATLAB 调用函数程序function y1=fl(x)y1=15/(2*x+1);② 将MATLAB 主程序和调用函数程序分别命名liti112.m 和fl.m ，分别保存为M 文件，然后在MATLAB 工作窗口输入命令>> [k,juecha,xiangcha,xk]= liti112(2,2.5,100) ③运行后输出计算结果列入表1–3和表 1-4中.④将算法2的MATLAB 调用函数程序的函数分别用y1=15-2*x^2和y1=x-(2*x^2+x-15)/(4*x+1)代替，得到算法1和算法3的调用函数程序，将其保存，运行后将三种算法的前8个迭代值821,,,x x x 列在一起（见表 1-3），进行比较.将三种算法的821,,,x x x 对应的绝对误差和相对误差的值列在一起（见表 1-4），进行比较.1.4 数值计算中应注意的问题例1.4.1 求数）181(71915-+⨯=-x 的近似值. 解 （1）直接用MATLAB 命令>> x=(7^15)*(sqrt(1+8^(-19))-1)运行后输出结果x = 0问题出现在两个相近的数1981-+与1相减时，计算机运行程序>>sqrt(1+8^(-19))-1运行后输出结果ans = 0由于计算机硬件只支持有限位机器数的运算，因此在计算中可能引入和传播舍入误差.因为有效数字的严重损失，导致输出18119-+-的结果为0，计算机不能再与数157继续进行真实的计算，所以，最后输出的结果与x 的精确值不符.（2）如果化为18187)181(71919151915++⨯=-+⨯=---x ，再用MATLAB 命令>> x=(7^15)*( (8^(-19))/(sqrt(1+8^(-19))+1))运行后输出结果x = 1.6471e-005 这是因为18119-+-化为18181919++--后，计算机运行程序>> x= (8^(-19))/(sqrt(1+8^(-19))+1)运行后的结果为x =3.4694e-018 由于有效数字的损失甚少，所以运算的结果-18103.4694⨯再与157继续计算，最后输出的结果与x 的精确值相差无几.例1.4.2 求数)13030ln(2--=y 的近似值. 解 （1）直接用MATLAB 程序>> x=30;x1= sqrt(x^2-1)运行后输出结果x1 = 29.9833 输入MATLAB 程序>> x=30; x1=29.9833;y=log(x-x1)运行后输出结果y = -4.0923（2）因为)13030ln(2--中的30=x 很大，如果采用倒数变换法111221-+=--=x x x x z ，即130301ln)13030ln(22-+=--)190030ln(-+-=.输入MATLAB 程序>> x=30;y=-log(x+sqrt(x^2-1))运行后输出结果y = -4.0941（3）输入MA TLAB 程序>> x=30; y=log(x-sqrt(x^2-1))运行后输出结果y = -4.0941 可见，（2）计算的近似值比（1）的误差小.参加计算的数，有时数量级相差很大.如果不注意采取相应的措施，在它们的加减法运算中，绝对值很小的那个数经常会被绝对值较大的那个数“吃掉”，不能发挥其作用，造成计算结果失真.例1.4.4 请在16位十进制数值精度计算机上利用软件MATLAB 计算下面的两个数0.30.1111111111111111*++=x 和0.30.11111111111111111*++=y将计算结果与准确值比较，解释计算结果.解 在MATLAB 工作窗口输入下面程序>> x=111111*********+0.1+0.3, y=1111111111111111+0.1+0.3运行后输出结果x = 1.111111*********e+014，y =1.111111*********e+015 从输出的结果可以看出，x *x =，而y *y ≠.为什么*y 仅仅比*x 多一位1，而y *y ≠呢？这是因为计算机进行运算时，首先要把参加运算的数写成绝对值小于1而“阶码”相同的数，这一过程称为数的“对阶”.在16位十进制数值精度计算机上利用软件MATLAB 计算这两个数，把运算的数*x 写成浮点规格化形式为，151515*103000**********.0001010000000000000.000100111111111111111.0⨯+⨯+⨯=x在16位十进制数值精度计算机上，三项的数都表示为小数点后面16位数字的数与1510之积，所以，计算机没有对数进行截断，而是按原来的三个数进行计算.因此，计算的结果x *x =.而161616*10030000000000000.00010010000000000000.000101111111111111111.0⨯+⨯+⨯=y三项的数都表示写成绝对值小于1而“阶码”都为1610的数以后，第一项的纯小数的小数点后面有16位数字.但是，后两项数的纯小数的小数点后面有17位数字，超过了16位十进制数值精度计算机的存储量，计算机对后两项的数都进行截断最后一位，即后两项的数都是16位机上的零，再进行计算，所以计算结果与实际不符.。

第一章绪论误差来源：模型误差、观测误差、截断误差（方法误差）、舍入误差是的绝对误差,是的误差，为的绝对误差限(或误差限）为的相对误差，当较小时，令相对误差绝对值得上限称为相对误差限记为：即:绝对误差有量纲，而相对误差无量纲若近似值的绝对误差限为某一位上的半个单位，且该位直到的第一位非零数字共有n位，则称近似值有n位有效数字，或说精确到该位。

例:设x==3。

1415926…那么,则有效数字为1位，即个位上的3，或说精确到个位.科学计数法：记有n位有效数字，精确到。

由有效数字求相对误差限：设近似值有n位有效数字，则其相对误差限为由相对误差限求有效数字：设近似值的相对误差限为为则它有n位有效数字令1.x+y近似值为和的误差（限）等于误差(限）的和2.x-y近似值为3.xy近似值为4.1．避免两相近数相减2．避免用绝对值很小的数作除数3．避免大数吃小数4．尽量减少计算工作量第二章非线性方程求根1。

逐步搜索法设f (a) <0, f (b）〉 0，有根区间为 (a, b）,从x0=a出发，按某个预定步长（例如h=（b-a）/N）一步一步向右跨，每跨一步进行一次根的搜索，即判别f(x k）=f（a+kh）的符号，若f（x k)〉0(而f（x k-1）<0)，则有根区间缩小为［x k-1，x k] （若f（x k）=0,x k即为所求根)，然后从x k—1出发，把搜索步长再缩小，重复上面步骤，直到满足精度：｜x k—x k-1｜< 为止，此时取x*≈(x k+x k-1）/2作为近似根.2。

二分法设f（x）的有根区间为［a,b]= [a0，b0], f（a）<0，f（b)〉0。

将［a0,b0]对分,中点x0= （(a0+b0)/2），计算f（x0）。

3.比例法一般地,设 [a k，b k］为有根区间，过（a k，f（a k）)、 (b k, f（b k）)作直线，与x轴交于一点x k，则:1.试位法每次迭代比二分法多算一次乘法，而且不保证收敛.2。

第一章数值计算中的误差分析数值计算方法(也称计算方法，数值方法)：是研究科学与工程技术中数学问题的数值解及其理论的一个数学分支，它的涉及面很广，涉及代数、微积分、微分方程数值解等问题。

●数值计算方法的主要任务：研究适合于在计算机上使用的数值计算方法及与此相关的理论，如方法的收敛性、稳定性以及误差分析等，此外，还要根据计算机的特点研究计算时间最短、需要计算机内存最少等计算方法问题．●数值计算主要过程：实际问题→建立数学模型→设计高效、可靠的数值计算方法→程序设计→上机计算求出结果。

●数值计算方法不同于纯数学：它既具有数学的抽象性与严格性，又具有应用的广泛性与实际试验的技术性，它是一门与计算机紧密结合的实用性很强的有着自身研究方法与理论系统的计算数学课程。

●数值计算方法的特点：应提供能让计算机直接处理的，包括加减乘除运算和逻辑运算及具有完整解题步骤的，切实可行的有效算法与程序，它可用框图、算法语言、数学语言或自然语言来描述，并有可靠的理论分析，能逼近且达到精度要求，对近似算法应保证收敛性和数值稳定性、进行必要的误差分析。

此外，还要注意算法能否在计算机上实现，应避免因数值方法选用不当、程序设计不合理而导致超过计算机的存贮能力，或导致计算结果精度不高等．根据“数值计算”的特点，首先应注意掌握数值计算方法的基本原理和思想，注意方法处理的技巧及其与计算机的密切结合，重视误差分析、收敛性及稳定性的基本理论；其次还要注意方法的使用条件，通过各种方法的比较，了解各种方法的异同及优缺点。

§1.1 误差的来源在数值计算过程中，估计计算结果的精确度是十分重要的工作，而影响精确度的因素是各种各样的误差，它们可分为两大类：一类称为“过失误差”，它一般是由人为造成的，这是可以避免的，故在数值计算中我们不讨论它；而另一类称为“非过失误差”，这在“数值计算”中往往是无法避免的，也是我们要研究的。

按照它们的来源，误差可分为以下四种：模型误差、观测误差、截断误差、舍入误差。