计算物理2013秋3幻灯片.pptVIP

1. 线性代数方程组求解 1.1 利用左除运算符的直接解法 对于线性方程组Ax=b，可以利用左除运算符“\”求解： x=A\b 例3-1 用直接解法求解下列线性方程组。 命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]; x=A\b 1.2 利用X=inv(A)*b语句的直接解法 方程AX＝b可改写为 X＝A－1 b 因此可用 X=inv(A)*b 解出。 例3-2 A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]; X= X=inv(A)*b 1.3最简行阶梯法解线性代数方程组 问题的提出 U=rref([A,b]) 程序运行的结果为： X1 = 15.58 X2 = 14.70 X3 = -8.20 X4 = 8.66 为什么要提出这种新的计算方法？ 把上例中第四个方程改为： ，求其解。 解：输入新参数 A=[6,1,6,-6;1,-1,9,9;-2,4,0,-4;4,2,7,-778/222]; b=[7;5;-7;877/222]; 键入U=rref([A,b])，得到 这个最简行阶梯形式说明原来的方程组是欠定的。 欠定方程组解的特点 它等价于下列方程组： 这是一组包括四个变量的三个有效方程。因此没有唯一的解。其中x4可以任意设定，即可以看作任意常数c，： 代入不同的c可以得到不同的解，因此欠定方程组有无数个解。这些解组成一根在空间中的直线。 用前两种方法试试？ 不相容方程组的示例 如果把第四个方程的右端常数仍取为 -9，则其行阶梯变换的结果为： 最后一个方程成了一个矛盾方程0=1。这说明方程组不相容，无解。 由此也可以看出，线性方程组求解最好还是用行阶梯简化的方法。因为它可以给出线性方程组的特征，避免计算的盲目性。 2.线性方程组应用举例 2.1多项式插值问题 给出平面上4个点的坐标值如右表， (1)。求对它进行插值的三次多项式， (2)。求t=1.5处f的近似值。 (3)。如果要求此多项式多通过一点(-1,5)，求其系数。 解：用多项式 来插值，令它在四点上的值与表中相同，，得到 多项式插值问题的解(1) 这个矩阵方程达到四阶，应该用计算机辅助求解了，编出 MATLAB程序exn554 C=[1,0,0,0;1,1,1,1;1,2,4,8;1,3,9,27], b=[3;0;-1;6],U=rref([C,b]) 得到， 多项式为 多项式插值问题的解(2) （2）把t1=1. 5代入此多项式，键入： f1=3-2*1.5-2*1.5^2+1.5^3 得到：f1=-1.125。 可以用以下语句画出插值图形， ezplot(t^3-2*t^2-2*t+3,[-1,4]), grid on,hold on plot(0:3,[3,0,-1,6],*) plot(1.5,-1.125,or) 可以得到图5-39，曲线通过了图中四个给定的插值点（用*号表示），圆圈为f(1.5)的位置。 多项式插值问题的解(3) (3)。若要此三次多项式多通过一点(-1,5)，将此点坐标代入后得，方程组就成为五个方程四个未知数，很可能是矛盾的，要靠计算来判断。用以下程序： A1=[1,0,0,0;1,1,1,1;1,2,4,8;1,3,9,27;1,-1,1,-1], b1=[3;0;-1;6;5],U01=rref([A1,b1]) 得到 多项式插值问题的解(4) 注意到系数增广矩阵最后一列是常数项，得出右方的同解方程组。其最后一个方程成为0=1,说明方程组是不相容的，属于超定方程，无通常意义下的解。 用MATLAB的pinv函数可以求出它的最小二乘解，键入： a1=pinv(A1)*b1，得： （pinv是一个求广义逆的函数） 多项式成为 多项式插值问题的解(5) 画出它的图形，可以继续键入 t=-1:0.1:4;plot(t,a1(1)+a1(2)*t+a1(3)*t.^2+a1(4)*t.^3,:r) plot(-1,5,x) 此曲线也用虚线画在图上。它不通过给定的5个点中的任何一个，说明都有误差。误差向量E及五项误差的平方和EE可以用矩阵方程计算如下 E=A1*a1-b1，EE=norm(E) 求得，结果为