文章目录

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

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

162页第12题（1）

用Newton法求 x3−x2−x−1=0x^3 - x^2 - x - 1 = 0x3−x2−x−1=0 的根，取x0=2x_0 = 2x0​=2，要求xk−xk−1<10−5x_k - x_{k - 1} < 10^{-5}xk​−xk−1​<10−5。

解：Newton迭代法的公式为xk+1=xk−f(xk)f′(xk),(k=0,1,2,...)x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}, (k = 0, 1, 2, ...)xk+1​=xk​−f′(xk​)f(xk​)​,(k=0,1,2,...)，

代码如下：

from sympy import *def g1(x):return x ** 3 - x ** 2 - x - 1def 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 = 1while (not abs(y - x) < 10 ** -5):print("第" + str(i) + "次迭代：当前解：" + str(y) + '\t' + "上一步解：" + str(x) + "\t" + "误差：" + str(abs(y - x)))x = yy = f(x)i = i + 1print("第" + 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)=x3+2x2+10x−20=0,f(x) = x^3 + 2x^2 + 10x - 20 = 0, f(x)=x3+2x2+10x−20=0,
并得出一个根α=1.36880817\alpha = 1.36880817α=1.36880817，但当时无人知道他用了什么方法，这个结果在当时是个非常著名的结果，请你构造一种简单迭代来验证此结果。

解：经计算：
f(1.3)=−1.423,f(1.4)=0.664,f(1.3)⋅f(1.4)<0,当1.3≤x≤1.4时，f′(x)>0,即f(x)在[1.3,1.4]上有且仅有一个解，则使用Newton法进行求证：f(1.3) = -1.423, f(1.4) = 0.664, f(1.3)·f(1.4) < 0, 当1.3 \leq x \leq 1.4时， f'(x) > 0,即f(x)在[1.3, 1.4]上有且仅有一个解，则使用Newton法进行求证： f(1.3)=−1.423,f(1.4)=0.664,f(1.3)⋅f(1.4)<0,当1.3≤x≤1.4时，f′(x)>0,即f(x)在[1.3,1.4]上有且仅有一个解，则使用Newton法进行求证：
代码如下：

from sympy import *def g1(x):return x ** 3 + 2 * x ** 2 + 10 * x - 20def 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 = 1while (not abs(y - x) < 10 ** -8):print("第" + str(i) + "次迭代：当前解：" + str(y) + '\t' + "上一步解：" + str(x) + "\t" + "误差：" + str(abs(y - x)))x = yy = f(x)i = i + 1print("第" + str(i) + "次迭代：当前解：" + str(y) + '\t' + "上一步解：" + str(x) + "\t" + "误差：" + str(abs(y - x)))if __name__ == '__main__':solve(1.3)

运行结果如图所示：

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

216页第12题

（数值实验题）人造地球卫星的轨道可视为平面上的椭圆，地心位于椭圆的一个焦点处。已知一颗人造地球卫星近地点距地球表面439km，远地点距地球表面2384km，地球半径为6371km。求该卫星的轨道长度。
解：长半轴a=(2384+439+6371∗2)/2=7782.5kma = (2384 + 439 + 6371 * 2) / 2 = 7782.5kma=(2384+439+6371∗2)/2=7782.5km，半个焦距c=a−439−6371=972.5kmc = a - 439 - 6371 = 972.5kmc=a−439−6371=972.5km，则短半轴b=a2−c2=7721.5kmb = \sqrt{a^2 - c^2} = 7721.5kmb=a2−c2​=7721.5km。设参数方程
{x=7782.5∗cos⁡ty=7721.5∗sin⁡t\left\{ \begin{array}{l} x = 7782.5 * \cos {t}\\ y = 7721.5 * \sin {t} \end{array} \right. {x=7782.5∗costy=7721.5∗sint​
其中t∈[0,2π]t \in [0, 2\pi]t∈[0,2π]。根据弧长公式：S=∫αβx′2(t)+y′2(t)dtS = \int_{\alpha}^{\beta}\sqrt{x^{'2}(t) + y^{'2}(t)}dtS=∫αβ​x′2(t)+y′2(t)​dt可知，周长积分公式如下：
L=4∗∫0π2(7782.5∗sin⁡x)2+(7721.5∗cos⁡x)2dxL = 4 * \int_{0}^{\frac{\pi}{2}} \sqrt{(7782.5*\sin{x})^2 + (7721.5*\cos{x})^2} dx L=4∗∫02π​​(7782.5∗sinx)2+(7721.5∗cosx)2​dx
使用GaussGaussGauss型积分公式计算上述积分，代码如下：

from sympy import *
import numpy as npdef f(x):return sqrt(((7782.5 * sin(x)) ** 2) + ((7721.5 * cos(x)) ** 2))def formula_cal():a = 7782.5b = 7721.5res = 2 * 3.1415926 * b + 4 * (a - b)print("公式计算轨道周长为：" + str(res) + "千米")def cal_integrate(i):# 权函数为1，i是需要增加的x的次数x = symbols('x')res = integrate(x ** i, (x, 0, 3.1415926 / 2))return resdef gauss_points(n):# 点的个数减一A = Matrix(np.zeros((n + 2, n + 2)))for i in range(0, n + 2):for j in range(0, n + 1):A[i, j] = cal_integrate(i + j)x = symbols('x')for i in range(0, n + 2):A[i, n + 1] = x ** iresult = solve(det(A), x)return resultdef gauss_cosfficient(n):# 点的个数减一b = gauss_points(n)# 拿到n个求积节点A = []x1 = []for i in range(0, len(b)):x1.append(f(b[i]))x = symbols('x')expr = 1.0 # 权函数for j in range(0, len(b)):if(j != i):expr = expr * ((x - b[j]) ** 2) / ((b[i] - b[j]) ** 2)A.append(integrate(expr, (x, 0, 3.1415926 / 2)))gauss_solve(x1, A)def gauss_solve(x1, A):n = len(x1)result = 0.0for i in range(0, n):result = result + x1[i] * A[i]print(str(n) + "点Gauss型积分求轨道周长为：" + str(result * 4) + "千米")if __name__ == '__main__':gauss_cosfficient(1)gauss_cosfficient(4)formula_cal()

运行结果如图所示：

其中计算轨道周长的公式为L=2πb+4(a−b)L = 2 \pi b + 4(a-b)L=2πb+4(a−b)，通过结果可以看出使用GaussGaussGauss型积分求解与使用公式求解的结果大致相等，正确性得到了验证。

《数值分析方法与应用》

一、基础知识部分

1、

设 SN=∑j=2N1j2−1S_N = \sum_{j = 2}^{N} \frac{1}{j^2 - 1}SN​=∑j=2N​j2−11​，其精确值为12(32−1N−1N+1)\frac{1}{2}(\frac{3}{2} - \frac{1}{N} - \frac{1}{N + 1})21​(23​−N1​−N+11​)。

（1）编制按从大到小的顺序 SN=122−1+132−1+⋅⋅⋅+1N2−1S_N = \frac{1}{2^2 - 1} + \frac{1}{3^2 - 1} + ··· + \frac{1}{N^2 - 1}SN​=22−11​+32−11​+⋅⋅⋅+N2−11​，计算SNS_NSN​的通用程序。

（2）编制按从小到大的顺序 SN=1N2−1+1(N−1)2−1+⋅⋅⋅+122−1S_N = \frac{1}{N^2 - 1} + \frac{1}{(N - 1)^2 - 1} + ··· + \frac{1}{2^2 - 1}SN​=N2−11​+(N−1)2−11​+⋅⋅⋅+22−11​，计算SNS_NSN​的通用程序。

（3）按两种顺序分别计算S102S_{10^2}S102​，S104S_{10^4}S104​，S106S_{10^6}S106​，并指出有效位数（编制程序时用单精度）。

（4）通过本上机题，你明白了什么。

解：代码如下：

import numpy as npdef f1(n):#从大到小res = np.float32(0.0)for i in range(2, n + 1):res = res + np.float32(1) / np.float32(i ** 2 - 1)print("N=" + str(n) + "，从大到小计算结果为" + str(np.float32(res)))def f2(n):#从小到大res = np.float32(0)for i in reversed(range(2, n + 1)):res = res + np.float32(1) / np.float32(i ** 2 - 1)print("N=" + str(n) + "，从小到大计算结果为" + str(np.float32(res)))def f3(n):#精确结果res = np.float32(3 * n * n - n - 2) / np.float32(4 * n * n + 4 * n)print("N=" + str(n) + "，精确结果为" + str(np.float32(res)))if __name__ == '__main__':f1(100)f2(100)f3(100)f1(10000)f2(10000)f3(10000)f1(1000000)f2(1000000)f3(1000000)

由于pythonpythonpython默认的都是双精度的浮点数，所以在这里我们用numpy.float32()来进行强制转换，将数据变成单精度的浮点数。

运行结果如图所示：

有效位数如下表所示：

从大到小 从小到大
10210^2102 7位 7位
10410^4104 4位 4位
10610^6106 3位 8位

通过上述实验数据可以看出：

（1）计算机存在舍入误差，有时会导致计算结果与精确结果略有出入；

（2）随着N的增大，有些分母会变得非常大，此时舍入误差的现象会变得很明显；

（3）此次算法使用从小到大的顺序进行得到的数据相对而言更精确，可以得到这样的启示：在数值计算时，要先分析不同算法对结果的影响，避免**“大数吃小数”**的现象，找出能得到更精确的结果的算法。

5、

首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用程序，程序包括输入多项式的系数以及给定点，输出函数值，利用编制的程序计算
p(x)=(x−2)9=x9−18x8+144x7−672x6+2016x5−4032x4+5376x3−4608x2+2304x−512p(x) = (x - 2)^9 = x^9 - 18x^8 + 144x^7 - 672x^6 + 2016x^5 - 4032x^4 + 5376x^3 -4608x^2 + 2304x - 512 p(x)=(x−2)9=x9−18x8+144x7−672x6+2016x5−4032x4+5376x3−4608x2+2304x−512
在x=2x = 2x=2邻域附近的值，画出p(x)p(x)p(x)在x∈[1.95,2.05]x \in [1.95, 2.05]x∈[1.95,2.05]上的图像。

解：代码如下：

import numpy as np
import matplotlib.pyplot as plt
def QJS(a, x):''':param a: 系数向量:param x: 给定点:return: 无'''res = 0.0if len(a) == 0:print("没有参数")else:res = a[0]for i in range(1, len(a)):res = res * x + a[i]print("函数在x = " + str(x) + "时的值为：" + str(res))def graph_function(a):''':param a: 系数向量:return: 无'''x = np.arange(1.95, 2.05, 0.01)y = a[0] * x ** 9 + a[1] * x ** 8 + a[2] * x ** 7 + a[3] * x ** 6 + a[4] * x ** 5 + a[5] * x ** 4 + a[6] * x ** 3 + a[7] * x ** 2 + a[8] * x + a[9]plt.figure()plt.plot(x, y, linestyle = '-', color = 'red')plt.xlabel('x')plt.ylabel('y')ax = plt.gca()ax.spines['right'].set_color('none')ax.spines['top'].set_color('none')ax.xaxis.set_ticks_position('bottom')ax.yaxis.set_ticks_position('left')plt.show()if __name__ == '__main__':a = [1, -18, 144, -672, 2016, -4032, 5376, -4608, 2304, -512]# 多项式系数向量x_0 = 2.03# 给定点QJS(a, x_0)# 计算x = 2邻域附近的值graph_function(a)# p(x)在[1.95, 2.05]上的图像

计算x=2x = 2x=2邻域附近的值（此处取2.03）：

画出p(x)p(x)p(x)在x∈[1.95,2.05]x \in [1.95, 2.05]x∈[1.95,2.05]上的图像：

二、线性方程组求解

2、

编制程序求解矩阵A的Cholesky分解，并用程序求解方程组Ax=bAx = bAx=b，其中
A=[71−511927−527−117−19],b=(13,−9,6,0)TA = \left[\begin{array}{r} 7 & 1 & -5 & 1\\ 1 & 9 & 2 & 7\\ -5 & 2 & 7 & -1\\ 1 & 7 & -1 & 9\\ \end{array}\right], b = (13, -9, 6, 0)^T A=⎣⎢⎢⎡​71−51​1927​−527−1​17−19​⎦⎥⎥⎤​,b=(13,−9,6,0)T

解：代码如下：

import numpy as np
from scipy import linalgdef cholesky_resolve(A):# Cholesky分解，返回L矩阵和L转置L = linalg.cholesky(A, lower = True)L_T = linalg.cholesky(A)return [L, L_T]def solve_lower(A, b):# 计算Ax = b，其中A是下三角矩阵，返回xn = np.shape(A)[0]x = np.zeros((n, 1))for i in range(0, n):res = 0.0for j in range(0, i):res = res + x[j] * A[i][j]x[i] = (b[i] - res) / A[i][i]return xdef solve_upper(A, b):# 计算Ax = b，其中A是上三角矩阵，返回xn = np.shape(A)[0]x = np.zeros((n, 1))for i in reversed(range(0, n)):res = 0.0for j in range(i + 1, n):res = res + x[j] * A[i][j]x[i] = (b[i] - res) / A[i][i]return xif __name__ == '__main__':A = np.array([[7, 1, -5, 1],[1, 9, 2, 7],[-5, 2, 7, -1],[1, 7, -1, 9]])b = np.array([13, -9, 6, 0])res = cholesky_resolve(A)L = res[0]L_T = res[1]y = solve_lower(L, b)x = solve_upper(L_T, y)print("x = \n", x)

运行结果如图所示：

本题运用Cholesky分解将对称正定矩阵A分解为A=LLTA = LL^TA=LLT，其中L是对角元均为正数的下三角矩阵，然后利用公式：
{Ly=bLTx=y\left\{ \begin{array}{l} Ly = b\\ L^Tx = y \end{array} \right. {Ly=bLTx=y​
进行求解。

6、

令HHH表示n×nn \times nn×n的Hilbert矩阵，其中(i,j)(i, j)(i,j)元素是1/(i+j−1)1/(i + j - 1)1/(i+j−1)，bbb是元素全为1的向量，用Gauss消去法求解Hx=bHx = bHx=b，其中取(1)n=2;(2)n=5;(3)n=10(1) n = 2;(2) n = 5;(3) n = 10(1)n=2;(2)n=5;(3)n=10。

解：代码如下：

import numpy as npdef create_augment(n):# 创建系数矩阵为Hilbert矩阵，b均为1的增广矩阵H = np.zeros((n, n + 1))for i in range(n):for j in range(n):H[i][j] = 1.0 / (i + j + 1)H[i][n] = 1.0return Hdef solve_upper(A, b):# 计算Ax = b，其中A是上三角矩阵，返回xn = np.shape(A)[0]x = np.zeros((n, 1))for i in reversed(range(0, n)):res = 0.0for j in range(i + 1, n):res = res + x[j] * A[i][j]x[i] = (b[i] - res) / A[i][i]return xdef gauss_solve(n):# 使用Gauss消元法将增广矩阵化成行阶梯型，并求解，返回xA = create_augment(n)for i in range(0, n - 1):for j in range(i + 1, n):alpha = A[j][i]/A[i][i]for k in range(i, n + 1):A[j][k] = A[j][k] - alpha * A[i][k]tempA = np.zeros((n, n))b = np.zeros((n, 1))for i in range(0, n):for j in range (0, n):tempA[i][j] = A[i][j]for i in range(0, n):b[i] = A[i][n]return solve_upper(tempA, b)if __name__ == '__main__':x_2 = gauss_solve(2)x_5 = gauss_solve(5)x_10 = gauss_solve(10)print("x_2 = \n", x_2)print("x_5 = \n", x_5)print("x_10 = \n", x_10)

运行结果如图所示：

x_2,x_5,x_10x\_2, x\_5, x\_10x_2,x_5,x_10分别是以2阶、5阶、10阶HilbertHilbertHilbert矩阵为系数矩阵，bbb全为1的线性方程组的解。

三、非线性方程组求解

1、

分别应用Newton迭代法和割线法计算（1）非线性方程2x3−5x+1=02x^3 - 5x + 1 = 02x3−5x+1=0在[1，2][1，2][1，2]上的一个根；（2）exsin⁡x=0e^x\sin x = 0exsinx=0在[−4,−3][-4, -3][−4,−3]上的一个根。

使用割线法时，由于迭代计算新的值需要前两个值，所以除了初始值之外，我们需要手动计算一个值；也可以先采用NewtonNewtonNewton迭代法先计算出来第二个值，然后代入到割线法中。

（1）使用NewtonNewtonNewton迭代法和割线法计算非线性方程2x3−5x+1=02x^3 - 5x + 1 = 02x3−5x+1=0在[1，2][1，2][1，2]上的一个根：

解：代码如下：

from sympy import *def g1(x):return 2.0 * x ** 3 - 5.0 * x + 1.0def newton(i):x = symbols("x")  # 符号x，自变量return i - g1(i) / diff(g1(x), x).subs('x', i)def secant(x_1, x_2):return x_2 - (g1(x_2) / (g1(x_2) - g1(x_1))) * (x_2 - x_1)def newton_solve(x):y = float(newton(x))i = 1while (not abs(y - x) < 10 ** -5):print("第" + str(i) + "次迭代：当前解：" + str(y) + '\t' + "上一步解：" + str(x) + "\t" + "误差：" + str(abs(y - x)))x = yy = newton(x)i = i + 1print("第" + str(i) + "次迭代：当前解：" + str(y) + '\t' + "上一步解：" + str(x) + "\t" + "误差：" + str(abs(y - x)))def secant_solve(x_1, x_2):a = x_1b = x_2i = 1# 第一步笔算，先走一步得到两个值while(not abs(b - a) < 10 ** -5):print("第" + str(i) + "次迭代：当前解：" + str(b) + '\t' + "上一步解：" + str(a) + "\t" + "误差：" + str(abs(b - a)))c = secant(a, b)a = bb = ci = i + 1print("第" + str(i) + "次迭代：当前解：" + str(b) + '\t' + "上一步解：" + str(a) + "\t" + "误差：" + str(abs(b - a)))if __name__ == '__main__':newton_solve(2.0)# 初始值选择2.0print("===================上面是Newton迭代法求解，下面是割线法求解=====================")secant_solve(2.0, 1.63157895)# 初始值选择2.0，第二个值是第一步采用切线法算出来的

运行结果如图所示：

（2）使用NewtonNewtonNewton迭代法和割线法计算非线性方程exsin⁡x=0e^x\sin x = 0exsinx=0在[−4,−3][-4, -3][−4,−3]上的一个根；

解：代码如下：

from sympy import *def g1(x):return E ** x * sin(x)def newton(i):x = symbols("x")  # 符号x，自变量return i - g1(i) / diff(g1(x), x).subs('x', i)def secant(x_1, x_2):return x_2 - (g1(x_2) / (g1(x_2) - g1(x_1))) * (x_2 - x_1)def newton_solve(x):y = float(newton(x))i = 1while (not abs(y - x) < 10 ** -5):print("第" + str(i) + "次迭代：当前解：" + str(y) + '\t' + "上一步解：" + str(x) + "\t" + "误差：" + str(abs(y - x)))x = yy = newton(x)i = i + 1print("第" + str(i) + "次迭代：当前解：" + str(y) + '\t' + "上一步解：" + str(x) + "\t" + "误差：" + str(abs(y - x)))def secant_solve(x_1, x_2):a = x_1b = x_2i = 1# 第一步笔算，先走一步得到两个值while(not abs(b - a) < 10 ** -5):print("第" + str(i) + "次迭代：当前解：" + str(b) + '\t' + "上一步解：" + str(a) + "\t" + "误差：" + str(abs(b - a)))c = secant(a, b)a = bb = ci = i + 1print("第" + str(i) + "次迭代：当前解：" + str(b) + '\t' + "上一步解：" + str(a) + "\t" + "误差：" + str(abs(b - a)))if __name__ == '__main__':newton_solve(-3.0)# 初始值选择-3.0print("===================上面是Newton迭代法求解，下面是割线法求解=====================")secant_solve(-3.0, -3.12476213)# 初始值选择-3.0，第二个值是第一步采用切线法算出来的

运行结果如图所示：

4、

用NewtonNewtonNewton迭代法求解方程x3+x2+x−3=0x^3 + x^2 + x - 3 = 0x3+x2+x−3=0的根，初值选择x0=−0.7x_0 = -0.7x0​=−0.7，迭代7步并与真值x∗=1x^* = 1x∗=1相比较，并列出数据表（表1）。

数据表(表 1)

iii xix_ixi​ ei=∣xi−x∗∣e_i = \vert x_i - x^* \vertei​=∣xi​−x∗∣ eiei−12\frac{e_i}{e_{i - 1}^2}ei−12​ei​​
0 -0.7 0.3 无
1 2.6205607476635517 1.6205607476635517 0.560747663551402
2 1.70844019026256 0.708440190262561 0.269756898741237
3 1.20637869801811 0.206378698018109 0.411205094190995
4 1.02416166412510 0.0241616641251030 0.567279521785564
5 1.00038149112102 0.000381491121020261 0.653477665330685
6 1.00000009699281 9.69928142247056e-8 0.666454786687556
7 1.00000000000001 6.21724893790088e-15 0.660874714617131

解：代码如下：

from sympy import *def g1(x):return x ** 3 + x ** 2 + x - 3.0def newton(i):x = symbols("x")  # 符号x，自变量return i - g1(i) / diff(g1(x), x).subs('x', i)def newton_solve(x):y = float(newton(x))i = 1while (i < 7):print("第" + str(i) + "次迭代：当前解：" + str(y) + '\t' + "误差：" + str(abs(y - 1.0)) + "\t" + "平方收敛系数" + str(abs(y - 1.0) / (abs(x - 1.0) * abs(x - 1.0))))x = yy = newton(x)i = i + 1print("第" + str(i) + "次迭代：当前解：" + str(y) + '\t' + "误差：" + str(abs(y - 1.0)) + "\t" + "平方收敛系数" + str(abs(y - 1.0) / (abs(x - 1.0) * abs(x - 1.0))))if __name__ == '__main__':newton_solve(-0.7)

运行结果如图所示：

四、插值与逼近

1、

已知函数f(x)=11+x2f(x) = \frac{1}{1 + x^2}f(x)=1+x21​，在[−5,5][-5, 5][−5,5]上分别取2,1,122, 1, \frac{1}{2}2,1,21​为单位长度的等距节点作为插值节点，用Lagrange方法插值，并把原函数图与插值函数图比较，观察插值效果。

解：代码如下：

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
from pylab import *def g1(x):return 1.0 / (x * x + 1.0)def create_table(x):# 以x为单位长度划分区间A = np.zeros((int(10 / x + 1), 2))i = 0flag = -5.0while (flag <= 5):A[i][0] = flagA[i][1] = g1(flag)i = i + 1flag = flag + xreturn Adef lagrange(x):# 求插值多项式A = create_table(x)m = shape(A)[0]x = symbols('x')y = 0.0 * xfor i in range(0, m):temp = x ** 0.0for j in range(0, m):if(j == i):continuetemp = temp * (x - A[j][0]) / (A[i][0] - A[j][0])y = y + A[i][1] * tempy = expand(y)return ydef draw_graph():'''由于四个函数图像放在一张图里面绘制出来的图像很不直观，所以笔者选择绘制三张图象,每幅图像里面都包含一个原函数和一个求出来的Lagrange插值多项式其中红色图像是原函数，蓝色图像是Lagrange插值多项式'''mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体（解决中文无法显示的问题）mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像时负号“-”显示方块的问题x = np.arange(-5.0, 5.0, 0.5)y_0 = 1.0 / (x * x + 1)y_1 = 0.00192307692307692 * x ** 4 - 0.0692307692307692 * x ** 2 - 5.55111512312578e-17 * x + 0.567307692307692y_2 = -2.26244343891403e-5 * x ** 10 - 1.70465380634928e-20 * x ** 9 + 0.00126696832579186 * x ** 8 - 9.14795583034644e-20 * x ** 7 - 0.0244117647058824 * x ** 6 + 1.90819582357449e-17 * x ** 5 + 0.19737556561086 * x ** 4 + 6.96057794735694e-17 * x ** 3 - 0.67420814479638 * x ** 2 - 1.52764086103208e-16 * x + 1.0y_3 = 2.72817068149799e-9 * x ** 20 - 4.54173855079897e-23 * x ** 19 - 2.65314598775676e-7 * x ** 18 - 6.38649607339942e-20 * x ** 17 + 1.07425130797328e-5 * x ** 16 + 5.45325105398239e-18 * x ** 15 - 0.000236412102809865 * x ** 14 - 1.08589623838001e-16 * x ** 13 + 0.00310184793200539 * x ** 12 - 4.45315713557687e-15 * x ** 11 - 0.0251135266738683 * x ** 10 - 1.78681939036474e-15 * x ** 9 + 0.126252909857137 * x ** 8 - 2.24466712578364e-14 * x ** 7 - 0.391630076762948 * x ** 6 + 4.11996825544492e-18 * x ** 5 + 0.753353962815142 * x ** 4 - 8.79743326798188e-15 * x ** 3 - 0.965739184991246 * x ** 2 - 7.99057001121817e-17 * x + 1.0plt.figure()plt.plot(x, y_0, linestyle='-', color='red', label='原函数')plt.plot(x, y_1, linestyle = '-', color = 'blue', label = '间隔为2')#plt.plot(x, y_2, linestyle = '-', color = 'blue', label = '间隔为1')#plt.plot(x, y_3, linestyle = '-', color = 'blue', label = '间隔为0.5')plt.xlabel('x')plt.ylabel('y')plt.legend(loc=0)ax = plt.gca()ax.spines['right'].set_color('none')ax.spines['top'].set_color('none')ax.xaxis.set_ticks_position('bottom')ax.yaxis.set_ticks_position('left')plt.show()if __name__ == '__main__':y_1 = lagrange(2)y_2 = lagrange(1)y_3 = lagrange(0.5)print("y_1 = ", y_1)print("y_2 = ", y_2)print("y_3 = ", y_3)draw_graph()

运行结果如图所示：

由于受到宽度限制，截图并没有显示出全部结果，下方采用文本复制的方式给出最终结果：

y_1 = 0.00192307692307692 * x ** 4 - 0.0692307692307692 * x ** 2 - 5.55111512312578e-17 * x + 0.567307692307692
y_2 = -2.26244343891403e-5 * x ** 10 - 1.70465380634928e-20 * x ** 9 + 0.00126696832579186 * x ** 8 - 9.14795583034644e-20 * x ** 7 - 0.0244117647058824 * x ** 6 + 1.90819582357449e-17 * x ** 5 + 0.19737556561086 * x ** 4 + 6.96057794735694e-17 * x ** 3 - 0.67420814479638 * x ** 2 - 1.52764086103208e-16 * x + 1.0
y_3 = 2.72817068149799e-9 * x ** 20 - 4.54173855079897e-23 * x ** 19 - 2.65314598775676e-7 * x ** 18 - 6.38649607339942e-20 * x ** 17 + 1.07425130797328e-5 * x ** 16 + 5.45325105398239e-18 * x ** 15 - 0.000236412102809865 * x ** 14 - 1.08589623838001e-16 * x ** 13 + 0.00310184793200539 * x ** 12 - 4.45315713557687e-15 * x ** 11 - 0.0251135266738683 * x ** 10 - 1.78681939036474e-15 * x ** 9 + 0.126252909857137 * x ** 8 - 2.24466712578364e-14 * x ** 7 - 0.391630076762948 * x ** 6 + 4.11996825544492e-18 * x ** 5 + 0.753353962815142 * x ** 4 - 8.79743326798188e-15 * x ** 3 - 0.965739184991246 * x ** 2 - 7.99057001121817e-17 * x + 1.0

绘制的函数图像如下图所示：


通过对比上面三幅函数图像可知，插值得到的多项式在已知点的函数值与原函数对应的函数值是严格相等的，则这就意味着，在一定范围内，这样的函数点越多，得到的图像就与原函数越逼近。

5、

考虑函数f(x)=sin⁡πx,x∈[0,1]f(x) = \sin \pi x, x \in [0, 1]f(x)=sinπx,x∈[0,1]。用等距节点作f(x)f(x)f(x)的NewtonNewtonNewton插值，画出插值多项式以及f(x)f(x)f(x)的图像，观察收敛性。

解：代码如下：

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
from pylab import *def g1(x):return sin(pi * x)def create_table(x):# 以x为单位长度划分区间A = np.zeros((int(1.0 / x + 1), int(1.0 / x + 1) + 1))i = 0flag = 0.0while (flag <= 1.0):A[i][0] = flagA[i][1] = g1(flag)i = i + 1flag = flag + xreturn Adef lagrange(x):# 求插值多项式A = create_table(x)m = shape(A)[0]n = shape(A)[1]# 求矩阵for i in range(2, n):for j in range(i - 1, m):A[j][i] = (A[j][i-1] - A[j-1][i-1]) / (A[j][0] - A[j-i+1][0])# 生成Newton多项式x = symbols('x')y = A[0][1] * x ** 0for i in range(1, m):temp = x ** 0.0for j in range(0, i):temp = temp * (x - A[j][0])y = y + A[i][i+1] * tempy = expand(y)return ydef draw_graph():mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体（解决中文无法显示的问题）mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像时负号“-”显示方块的问题x = np.arange(-0.1, 1.1, 0.1)y_0 = sin(pi * x)y_1 = 3.61347072168978 * x ** 4 - 7.22694144337956 * x ** 3 + 0.517968210332188 * x ** 2 + 3.09550251135759 * xplt.figure()plt.plot(x, y_0, linestyle='-', color='red', label='原函数')plt.plot(x, y_1, linestyle = '-', color = 'blue', label = 'Newton插值多项式')plt.xlabel('x')plt.ylabel('y')plt.legend(loc=0)ax = plt.gca()ax.spines['right'].set_color('none')ax.spines['top'].set_color('none')ax.xaxis.set_ticks_position('bottom')ax.yaxis.set_ticks_position('left')plt.show()if __name__ == '__main__':y = lagrange(0.2)# 节点间距为0.2print("y = ", y)draw_graph()

运行结果如图所示：

上方的函数即为NewtonNewtonNewton插值生成的多项式函数，当前节点间距为0.2，生成函数与原函数图像对比（红色曲线为原函数，蓝色曲线为Newton多项式）如图所示：

可见在[0,1]这个区间上，生成多项式的收敛性还是比较好的。下面附上几个不同间距的生成多项式与原函数对比：

节点间距为0.5

Newton插值多项式：y = -4.0*x**2 + 4.0*x

节点间距为0.1

Newton插值多项式：y =  -0.024766311852928*x**10 + 0.123831559266328*x**9 - 0.0434827562171139*x**8 - 0.569058330736058*x**7 - 0.0143385072780853*x**6 + 2.55481222833824*x**5 - 0.00100448543439632*x**4 - 5.16757591409148*x**3 - 1.04708062454788e-5*x**2 + 3.14159298881174*x

可见节点间距为0.5时收敛性一般，间距为0.2时收敛性有所好转，间距为0.1时收敛性比前两者都要好。

7、

编程计算三次样条SSS，满足S(0)=1,S(1)=3,S(2)=3,S(3)=4,S(4)=2S(0) = 1, S(1) = 3, S(2) = 3, S(3) = 4, S(4) = 2S(0)=1,S(1)=3,S(2)=3,S(3)=4,S(4)=2，其中边界条件S′′(0)=S′′(4)=0S''(0) = S''(4) = 0S′′(0)=S′′(4)=0。
解：使用第二类边界条件，代码如下：

import numpy as np
from sympy import *def generate_g0(xy):return (3.0 * (xy[1][1] - xy[1][0]) / (xy[0][1] - xy[0][0])) - ((xy[0][1] - xy[0][0]) / 2.0 * 0)def generate_g4(xy):return (3.0 * (xy[1][4] - xy[1][3]) / (xy[0][4] - xy[0][3])) - ((xy[0][4] - xy[0][3]) / 2.0 * 0)def spline_solve(xy):# 生成系数矩阵A = np.zeros((5,5))A[0][0] = 2.0A[4][4] = 2.0A[0][1] = 1.0A[4][3] = 1.0for i in range(1, 4):A[i][i - 1] = 0.5A[i][i] = 2A[i][i + 1] = 0.5# 生成列向量gg = np.zeros((5, 1))g[0][0] = generate_g0(xy)g[4][0] = generate_g4(xy)for i in range(1, 4):g[i][0] = 3.0 / 2.0 * (((xy[1][i + 1] - xy[1][i]) / (xy[0][i + 1] - xy[0][i])) + ((xy[1][i] - xy[1][i - 1]) / (xy[0][i] - xy[0][i - 1])))# 求解Am = gm = np.linalg.solve(A, g)# 生成多项式for i in range(0, 4):x = symbols('x')s = (1 + 2 * (x - xy[0][i])) * ((x - xy[0][i + 1]) ** 2) * xy[1][i] + (1 - 2 * (x - xy[0][i + 1])) * ((x - xy[0][i]) ** 2) * xy[1][i + 1] + (x - xy[0][i]) * ((x - xy[0][i + 1]) ** 2) * m[i][0] + (x - xy[0][i + 1]) * ((x - xy[0][i]) ** 2) * m[i + 1][0]s = expand(s)print("x属于[" + str(i) + ", " + str(i+1) + "]时，s(x) = ", s)if __name__ == '__main__':xy = np.array([[0.0, 1.0, 2.0, 3.0, 4.0], [1.0, 3.0, 3.0, 4.0, 2.0]])spline_solve(xy)

运行结果如图所示：

经验证，第二段函数和第一段函数之间的差大约等于1.9643∗(x−1)31.9643*(x - 1)^31.9643∗(x−1)3，第三段函数和第二段函数之间的差大约等于−2.8572∗(x−2)3-2.8572*(x - 2)^3−2.8572∗(x−2)3，第四段函数和第三段函数之间的差大约等于2.4643∗(x−3)32.4643*(x - 3)^32.4643∗(x−3)3，所以该分片函数是所求的三次样条函数。

五、数值积分

2、

用两点，三点和五点的Gauss型积分公式分别计算定积分，并与真值作比较。
S=∫0π2x2cos⁡xdxS = \int_{0}^{\frac{\pi}{2}} x^2\cos x dx S=∫02π​​x2cosxdx
解：为了便于观察结果，程序中所有的π\piπ均使用近似值3.14159263.14159263.1415926代替，否则结果中均带有π\piπ，不利于比较。

代码如下：

from sympy import *
import numpy as npdef f(x):return cos(x)def real_integrate():# 真值x = symbols('x')res = integrate((x ** 2) * cos(x), (x, 0, 3.1415926 / 2))print("真值为：" + str(res))def cal_integrate(i):# 权函数为x平方，i是还需要增加的次数x = symbols('x')res = integrate(x ** (2 + i), (x, 0, 3.1415926 / 2))return resdef gauss_points(n):# 点的个数减一A = Matrix(np.zeros((n + 2, n + 2)))for i in range(0, n + 2):for j in range(0, n + 1):A[i, j] = cal_integrate(i + j)x = symbols('x')for i in range(0, n + 2):A[i, n + 1] = x ** iresult = solve(det(A), x)return resultdef gauss_cosfficient(n):# 点的个数减一b = gauss_points(n)# 拿到n个求积节点A = []x1 = []for i in range(0, len(b)):x1.append(f(b[i]))x = symbols('x')expr = x ** 2for j in range(0, len(b)):if(j != i):expr = expr * ((x - b[j]) ** 2) / ((b[i] - b[j]) ** 2)A.append(integrate(expr, (x, 0, 3.1415926 / 2)))gauss_solve(x1, A)def gauss_solve(x1, A):n = len(x1)result = 0.0for i in range(0, n):result = result + x1[i] * A[i]print(str(n) + "点Gauss型积分为：" + str(result))if __name__ == '__main__':gauss_cosfficient(1)gauss_cosfficient(2)gauss_cosfficient(4)real_integrate()

运行结果如图所示：

其中3点GaussGaussGauss型积分得到的结果带有虚数，比较长，在下方另行给出：

3点Gauss型积分为：(0.518479879429716 - 1.13170196486494e-23*I)*(0.566819007009947 + 8.17802788116717e-24*I) + (0.116082467091191 + 5.82713089875272e-24*I)*(0.894546194955815 - 2.95783579853275e-24*I) + (0.114407705350449 + 1.31479879463766e-23*I)*(0.609026654797592 + 6.18118654397656e-23*I)

在此算法中，笔者采用的是按照HermiteHermiteHermite基函数得到的求积系数公式。从结果可以看出，有着最高代数精度的GaussGaussGauss型求积公式，在有两点的时候就已经很接近真实结果了。

六、微分方程数值解法

1、

已知常微分方程
{dudx=2xu+x2exx∈[1,2],u(1)=0\left\{ \begin{array}{l} \frac{du}{dx} = \frac{2}{x}u + x^2e^x\\ x \in [1, 2],u(1) = 0 \end{array} \right. {dxdu​=x2​u+x2exx∈[1,2],u(1)=0​
分别用Euler法，改进的Euler法，Runge-Kutta法去求解该方程，步长选为0.1,0.05,0.010.1, 0.05, 0.010.1,0.05,0.01。画图观察求解效果。
解：本题采用4级4阶经典Runge−KuttaRunge-KuttaRunge−Kutta法，代码如下：

import numpy as np
import math
import matplotlib.pyplot as pltdef u(x):return pow(x, 2) * (pow(math.e, x) - pow(math.e, 1))def du_dx(x, ux):return (2 / x) * ux + pow(x, 2) * pow(math.e, x)def real_value(a, b, h):t = np.arange(a + h, b + h, h)un = [0]for i in t:un.append(u(i))print("Real-Values:")print(un)return undef Euler(a, b, h):''':param a: 左边界:param b: 右边界:param h: 步长'''un = [0] # u(1)=0res = 0# 以步数划分区间，x=1除外 x=1时u=0t = np.arange(a, b + h, h)for i in range(len(t) - 1):res = res + h * du_dx(t[i], un[-1])un.append(res)print("Euler:")print(un)return undef improvedEuler(a, b, h):''':param a: 左边界:param b: 右边界:param h: 步长'''un = [0] # u(1)=0res = 0# 以步数划分区间，x=1除外 x=1时u=0t = np.arange(a, b + h, h)for i in range(len(t) - 1):tempRes = res + h * du_dx(t[i], un[-1])res = res + (h / 2) * (du_dx(t[i], un[-1]) + du_dx(t[i + 1], tempRes))un.append(res)print("Improved-Euler:")print(un)return undef runge_kutta(y, x, dx, f):'''y is unx is tndx is step_lengthf is du_dx'''k1 = f(x, y)k2 = f(x + 0.5 * dx, y + 0.5 * dx * k1)k3 = f(x + 0.5 * dx, y + 0.5 * dx * k2)k4 = f(x + dx, y + dx * k3)return y + dx * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0def cal_runge_kutta(a, b, h, y):''':param a: 左边界:param b: 右边界:param h: 步长:param y: u(1)=1'''ys, ts = [0], [1.0]while a <= b:y = runge_kutta(y, a, h, du_dx)a += hys.append(y)ts.append(a)print("Runge-Kutta:")print(ys)return ts, ysif __name__ == '__main__':h = 0.1# h = 0.05# h = 0.01# 存储真实结果e0 = real_value(1.0, 2.0, h)# 储存Euler结果e1 = Euler(1.0, 2.0, h)# 储存improvedEuler结果e2 = improvedEuler(1.0, 2.0, h)# 储存Runge-Kutta结果ts, ys = cal_runge_kutta(1.0, 2.0, h, 0)# 绘图横坐标plt.plot(ts, e0, label = 'Real-Value')plt.plot(ts, e1, label = 'Euler')plt.plot(ts, e2, label = 'Improved-Euler')plt.plot(ts, ys, label = 'Runge-Kutta')# 设置横坐标的文字说明plt.xlabel("value of x")# 设置纵坐标的文字说明plt.ylabel("value of u")plt.title('h = %s' % h)plt.legend()plt.show()

注：步长可通过修改主函数中的hhh进行修改。

运行结果如图所示：

上图选取的步长为0.1，得出的数据结果较为冗余，不在此处一一列出。通过观察图像可以得出更为简洁明确的结论：



通过观察上述三幅根据不同步长得出的图像，可以看出在步长较大的时候，具有一阶精度的EulerEulerEuler法得出的结果与真实值有较大出入，而改进的EulerEulerEuler法和Runge−KuttaRunge-KuttaRunge−Kutta法在步长较大的时候就与真实结果十分接近；随着步长的缩小，EulerEulerEuler法与真实值越来越逼近，而改进的EulerEulerEuler法和Runge−KuttaRunge-KuttaRunge−Kutta法始终都与真实值十分接近。

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

  1. 2022年秋季学期人工神经网络第四次作业

    说明: 本次作业是针对这学期经典神经网络中的内容,主要涵盖竞争神经网络课程内容相关的算法. 完成作业可以使用你所熟悉的编程语言和平台,比如 C,C++.MATLAB.Python等. 作业要求链接: ...

  2. 【20保研】大连理工大学软件学院2019年优秀大学生学术夏令营通知

    点击文末的阅读原文或者公众号界面左下角的保研夏令营或者公众号回复"夏令营"是计算机/软件等专业的所有保研夏令营信息集合,会一直更新的. 一.院系与学科介绍  大连理工大学软件工程学 ...

  3. 大连理工大学c语言第三次上机作业答案,大连理工大学软件学院C语言上机第五六章课后题...

    大连理工大学软件学院C语言上机第五六章课后题 五.1. #includeint main() { int a,b,c; float X,Y,Z; scanf("%d%d%d",&a ...

  4. 大连理工大学软件学院编译技术课程——MicroC词法分析上机实验

    大连理工大学软件学院编译技术课程--MicroC词法分析上机实验 题目 编写词法分析编译程序 实验目的:对循环语句和条件判断语句编写词法分析编译程序,只能通过一遍扫描完成. 实验要求: (1) 关键字 ...

  5. 大连理工大学计算机考研调剂,2020年大连理工大学软件学院考研调剂公告

    考研调剂公告各院校已经发布出来了,下面由出国留学网小编为你精心准备了"2020年大连理工大学软件学院考研调剂公告",持续关注本站将可以持续获取更多的考试资讯! 2020年大连理工大 ...

  6. 大连理工计算机组成实验,大连理工大学软件学院计算机组成原理实验报告

    <大连理工大学软件学院计算机组成原理实验报告>由会员分享,可在线阅读,更多相关<大连理工大学软件学院计算机组成原理实验报告(57页珍藏版)>请在人人文库网上搜索. 1.大连理工 ...

  7. 大连理工大学软件学院博客地址

    大连理工大学软件学院CSDN高校俱乐部博客地址: http://blog.csdn.net/u010455967 转载于:https://www.cnblogs.com/FlightButterfly ...

  8. 【调剂】985大连理工大学软件学院2020年硕士研究生调剂缺额与报名通知

    点击文末的阅读原文或者公众号界面左下角的调剂信息或者公众号回复"调剂"是计算机/软件等专业的所有调剂信息集合,会一直更新的. 根据<软件学院学院2020年硕士研究生复试录取办 ...

  9. 2018大连理工计算机考研分数线,2018年大连理工大学软件学院考研复试分数线

    根据<大连理工大学2018年全国硕士研究生招生复试办法>(大工研发[2018]3号)的精神,结合我院实际情况,本着公平.公正.公开的原则,坚持执行政策透明.规则公开.程序公正.结果公开.监 ...

  10. 大连理工计算机保研面试,保研经【大连理工大学软件学院吧】_百度贴吧

    该楼层疑似违规已被系统折叠 隐藏此楼查看此楼 还有一些学校计算机系是接受软件学院学生的,这个一定要问清楚,非常的重要.04届的时候,这个问题,一般教务员就能告诉你结果了,因为去年不要的,今年还是那样, ...

最新文章

  1. 手机应用软件测试的思路与要点
  2. 跟着石头哥哥学cocos2d-x（三）---2dx引擎中的内存管理模型
  3. CTF web题总结--爆破用户名密码
  4. 处理大并发之一 对异步非阻塞的理解
  5. 【Python】青少年蓝桥杯_每日一题_10.19_回文数和个数
  6. JavaScript实现ShellSort希尔排序算法（附完整源码）
  7. 学php要懂js吗,js要怎么学
  8. 小程序 mpvue input 文本控制
  9. 微信端自动授权登陆实现 - 无第三方库版
  10. Oracle之触发器
  11. 大数据智能学院的硕士论文怎么写_大数据智能营销笔记本怎么样
  12. iOS sign in with Apple 苹果ID登录
  13. hunter 网络空间测绘
  14. 解决unable to find valid certification path to requested target
  15. pip下载报错：pip._vendor.urllib3.exceptions.SSLError: [SSL: DECRYPTION_FAILED_OR_BAD_RECORD_MAC] decry
  16. 限速器校验合格范围_限速
  17. 学习笔记：Qt程序打包发布
  18. php微信支付na,PHP公众号支付宝支付实现
  19. 【Linux共享内存】
  20. 【体系结构系列】体系结构概述

热门文章

  1. GMAP.NET应用及离线地图加载
  2. 利用Python批量把flv文件转换成mp4文件
  3. 计算机应用的最广领域,从乡镇企业的从业人员数，我们可以看出：
  4. oracle快照点,Oracle快照（snapshot）管理
  5. 架设自己的邮件服务器
  6. canvas画出闪瞎眼的简单图形
  7. 汇编语言cf,of,sf,zf
  8. 使用空驱动消除设备管理器里面的未知设备
  9. eclipse下改变 匹配标签和匹配括号的颜色
  10. TS 中 as 用法