大连理工大学软件学院2022年秋季学期《矩阵与数值分析》上机作业

Posted 2022-12-23 谈瀛洲

文章目录

• 《计算机科学计算》第二版
• • 162页第12题（1）
  • 162页第16题
  • 216页第12题
• 《数值分析方法与应用》
• • 一、基础知识部分
  • • 1、
    • 5、
  • 二、线性方程组求解
  • • 2、
    • 6、
  • 三、非线性方程组求解
  • • 1、
    • 4、
  • 四、插值与逼近
  • • 1、
    • 5、
    • 7、
  • 五、数值积分
  • • 2、
  • 六、微分方程数值解法
  • • 1、

《计算机科学计算》第二版

162页第12题（1）

用Newton法求 $x^3-x^2-x-1=0$ 的根，取$x_0=2$，要求$x_k-x_{k-1}<10^{-5}$。

解：Newton迭代法的公式为$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)},(k=0,1,2,...)$，

代码如下：

from sympy import *

def g1(x):
    return x ** 3 - x ** 2 - x - 1

def f(i):
    x = symbols("x")  # 符号x，自变量
    return i - g1(i) / diff(g1(x), x).subs('x', i);

def solve(x):
    y = float(f(x))
    i = 1
    while (not abs(y - x) < 10 ** -5):
        print("第" + str(i) + "次迭代：当前解：" + str(y) + '\\t' + "上一步解：" + str(x) + "\\t" + "误差：" + str(abs(y - x)))
        x = y
        y = f(x)
        i = i + 1
    print("第" + str(i) + "次迭代：当前解：" + str(y) + '\\t' + "上一步解：" + str(x) + "\\t" + "误差：" + str(abs(y - x)))

if __name__ == '__main__':
    solve(2)

运行结果如图所示：

综上：方程的解为1.83928675521416

162页第16题

Heonardo于1225年研究了方程：
$$f(x)=x^3+2x^2+10x-20=0,$$
并得出一个根$\alpha =1.36880817$，但当时无人知道他用了什么方法，这个结果在当时是个非常著名的结果，请你构造一种简单迭代来验证此结果。

解：经计算：
$$f(1.3)=-1.423,f(1.4)=0.664,f(1.3)\cdot f(1.4)<0,当1.3\le x\le 1.4时，f^{\prime }(x)>0,即f(x)在[1.3,1.4]上有且仅有一个解，则使用Newton法进行求证：$$
代码如下：

from sympy import *

def g1(x):
    return x ** 3 + 2 * x ** 2 + 10 * x - 20

def f(i):
    x = symbols("x")  # 符号x，自变量
    return i - g1(i) / diff(g1(x), x).subs('x', i);

def solve(x):
    y = float(f(x))
    i = 1
    while (not abs(y - x) < 10 ** -8):
        print("第" + str(i) + "次迭代：当前解：" + str(y) + '\\t' + "上一步解：" + str(x) + "\\t" + "误差：" + str(abs(y - x)))
        x = y
        y = f(x)
        i = i + 1
    print("第" + str(i) + "次迭代：当前解：" + str(y) + '\\t' + "上一步解：" + str(x) + "\\t" + "误差：" + str(abs(y - x)))

if __name__ == '__main__':
    solve(1.3)

运行结果如图所示：

此处设置的相邻两次结果之差不大于$10^{-8}$，可见最后结果为：1.36880810782137，Heonardo的结果得以验证。

216页第12题

（数值实验题）人造地球卫星的轨道可视为平面上的椭圆，地心位于椭圆的一个焦点处。已知一颗人造地球卫星近地点距地球表面439km，远地点距地球表面2384km，地球半径为6371km。求该卫星的轨道长度。
解：长半轴$a=(2384+439+6371\ast 2)/2=7782.5km$，半个焦距$c=a-439-6371=972.5km$，则短半轴$b=\sqrt{a^2-c^2}=7721.5km$。设参数方程
$$\begin{array}{l}x=7782.5\ast \cos t\\ y=7721.5\ast \sin t\end{array}$$
其中$t\in [0,2\pi ]$。根据弧长公式：$S={\int }_{\alpha }^{\beta }\sqrt{x^{{}^{\prime }2}(t)+y^{{}^{\prime }2}(t)}dt$可知，周长积分公式如下：
L = 4 ∗ ∫ 0 π 2 ( 7782.5 ∗ sin ⁡ x ) 2 + ( 7721.5 ∗ cos ⁡ x ) 2 d x L = 4 * \\int_0^\\frac\\pi2 \\sqrt(7782.5*\\sinx)^2 + (7721.5