突然想起來可以做做數值分析的作業，于是把室友的數值分析作業拿過來練手，寫一篇博客分享一下，其實我這個菜鳥把程式寫復雜了，有很多可以簡化的地方，由于本菜鳥其它作業還沒寫完，就不去簡化了，大家可以自行改正啦，

數值分析上機題

首先說一下自己的疑惑，對于第一題python怎么實作對ln(x)直接呼叫求導呢？是直接用泰勒展開后對多項式求導嗎？第二題是用Newton迭代法，要求出迭代初始值在什么范圍可以得到收斂解，這里如果用迭代程式去實作的話也會比較麻煩，那有沒有更好的方法去求解呢？希望有大神可以幫忙留言解惑，多謝，

一、題1——對數的近似求解

1.題目描述

**題目：**這里偷下懶，直接貼出原題的截圖吧，

2.python實作

不多BB，程式在這：

import numpy as np
from sympy import * #用于求導積分等科學計算
import math as m

x = Symbol('x')#x 變數
t = Symbol('t')
y1 = 1/(1+x) #公式
y2 = -1/(1+x) #公式
y3 = 2/(1-x**2) #公式

def func(m):
    res = m
    for j in range(1,m):
        res *= j
    return res

def ln_Tyalor(y):
    Tl_expr = y * (t-x)
    for i in range(1, 10):
        fac = func(i+1)
        f_n = diff(y, x, i)
        Tl_expr += (f_n / fac)*(t-x)**(i+1)
    return Tl_expr.subs({x:0})

#print(ln_Tyalor(y1))
sexpr1 = ln_Tyalor(y1)
sexpr2 = ln_Tyalor(y2)
sexpr3 = ln_Tyalor(y3)
A = sexpr1.subs({t:1}).evalf()
B = sexpr2.subs({t:-1/2}).evalf()
C = sexpr3.subs({t:1/3}).evalf()
print('ln2的值：', m.log(2, m.e))
print('方程ln(1+x)的10階泰勒展開計算結果為：', A,'\n','估計誤差為：', abs(m.log(2, m.e)-A))
print('方程ln(1/(1+x))的10階泰勒展開計算結果為：', B,'\n','估計誤差為：', abs(m.log(2, m.e)-B))
print('方程ln((1+x)/(1-x))的10階泰勒展開計算結果為：', C,'\n','估計誤差為：', abs(m.log(2, m.e)-C))

3.輸出結果

ln2的值： 0.6931471805599453
方程ln(1+x)的10階泰勒展開計算結果為： 0.645634920634921 
 估計誤差為： 0.0475122599250246
方程ln(1/(1+x))的10階泰勒展開計算結果為： 0.693064856150793 
 估計誤差為： 8.23244091517905e-5
方程ln((1+x)/(1-x))的10階泰勒展開計算結果為： 0.693146047390827 
 估計誤差為： 1.13316911820593e-6

一、題2——方程求根的Newton法

1.題目描述

**題目：**還是截圖方便，

2.python實作

這個程式寫的不好，由于寫完后，exp（）函式老是報錯說：整型資料不是exp物件的屬性，改了之后也沒實作自己的想法，那就這樣吧，沒時間啦，大家自己搞吧，，，，，

import numpy as np
from sympy import * #用于求導積分等科學計算
import math as m

def f(x):
    return 3*x**2 - m.exp(x) #該方程有3個根，初步估計在[-1,0],[0,1],[3,4]

def fdiff(x):
    return 6*x - m.exp(x)

def g(x):
    return 6*x - m.exp(x)#該方程有3個根，初步估計在[0,1],[2,3]

def gdiff(x):
    return 6 - m.exp(x)

a = float(input('請輸入計算區間的下界a(浮點型）: '))
b = float(input('請輸入計算區間的上界b(浮點型）: '))
c = float(input('請輸入迭代初始值(浮點型）: '))
n = input('請輸入要求解的函式，f代表f(x),g代表g(x): ')

if n =='f':

    if f(a) * f(b)> 0:
        print('在此區間內函式有多個根或者無根')
    elif f(a) * f(b) == 0:
        print('f(a) = ', '%f'%f(a))
        print('f(b) = ', '%f'%f(b))
    else:
        fcount = 0
        y = c - f(c) / fdiff(c)
        while (abs(c - y) >= 0.5e-9) & (fdiff(c) != 0):
            x2 = c - f(c) / fdiff(c)
            y = c
            c = x2
            fcount += 1
        print('函式給出的求根區間是：', [a, b])
        print("函式f(x）的Newton迭代次數：%f,函式f(x）的迭代計算所得的根為:%.8f"%(fcount,c))

elif n =='g':

    if g(a) * g(b)> 0:
        print('在此區間內函式有多個根或者無根')
    elif g(a) * g(b) == 0:
        print('g(a) = ', '%f'%g(a))
        print('g(b) = ', '%f'%g(b))
    else:
        gcount = 0
        y = c - g(c) / gdiff(c)
        while (abs(c - y) >= 0.5e-9) & (gdiff(c) != 0):
            x2 = c - g(c) / gdiff(c)
            y = c
            c = x2
            gcount += 1
        print('函式給出的求根區間是：', [a, b])
        print("函式g(x）的Newton迭代次數：%f,函式g(x）的迭代計算所得的根為:%.8f"%(gcount,c))

3.輸出結果

這里說明一下，輸入的[a, b]是你想判斷該區間有沒有根；c是迭代初始值；f代表f（x)，g代表g(x)；這幾個引數請自己輸入，

請輸入計算區間的下界a(浮點型）: -1.0
請輸入計算區間的上界b(浮點型）: 4.0
請輸入迭代初始值(浮點型）: -3.0
請輸入要求解的函式，f代表f(x),g代表g(x): f
函式給出的求根區間是： [-1.0, 4.0]
函式f(x）的Newton迭代次數：7.000000,函式f(x）的迭代計算所得的根為:-0.45896227

一、題3——Hilbert矩陣的條件數

1.題目描述

**題目：**你懂的，

2.python實作

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['simhei']

n = int(input('請輸入Hilbert方陣的最大維數:' ))
def max_sum_rows(X):
    sum_row_list1 = []
    for i in range(len(X)):
        count = 0
        for j in range(len(X)):
            count += abs(X[i][j])
        sum_row_list1.append(count)
    return max(sum_row_list1)

inf_fanshu = []
Hilbert_cond = []
for i in range(1, n+1):
    X = 1./(np.arange(1, i+1) + np.arange(0, i)[:, np.newaxis])
    invX = np.linalg.inv(X)
    a1 = max_sum_rows(invX)
    a2 = max_sum_rows(X)
    inf_fanshu.append(a2)
    H_cond = a1 * a2
    Hilbert_cond.append(H_cond)

#計算10,20……100的無窮范數條件數
Hilbert_cond_test = []
for j in range(n):
    if (j+1)%10 == 0:
        Hilbert_cond_test.append(Hilbert_cond[j])
print('生成維數從10,20……100的Hilbert矩陣的行范數條件數：\n', Hilbert_cond_test)

plt.title('Hilbert矩陣維數與條件數對數之間的關系')
plt.plot((list(range(1,n+1))), np.log(Hilbert_cond),c ='b', marker='*',label='擬合曲線')
plt.legend()
plt.xlabel('矩陣維度n')
plt.xticks(np.arange(0, n+1, 5))
plt.yticks(np.arange(0, 55, 5))
plt.ylabel('log(cond)')
plt.show()

#求解Hilbert方程的解和對應的無窮條件數
r_A_A_acc_list = []
r_B_list = []
r_cond = []
r_B_cond = []
for i in range(1,n+1):
    A = np.ones((i,1))*1
    X = 1. / (np.arange(1, i + 1) + np.arange(0, i)[:, np.newaxis])
    B = X@A
    A_acc = np.linalg.inv(X)@B
    r_A_A_acc = A - A_acc
    r_B = B - X @ A_acc
    r_A_A_acc_list.append(r_A_A_acc)
    r_B_list.append(r_B)
    r_cond.append(abs(r_A_A_acc[:]).max()) #x-x_acc的無窮范數
    r_B_cond.append(abs(r_B[:]).max())#b-Hx_acc的無窮范數

print('維數為10,50,100時的x-x_acc計算結果：\n', r_A_A_acc_list[9] ,r_A_A_acc_list[49] , r_A_A_acc_list[99])
print('維數為10,50,100時的b-Hx_acc計算結果：\n',r_B_list[9], r_B_list[49], r_B_list[99])
print('維數為10,50,100時的x-x_acc矩陣的無窮條件數計算結果：\n', r_cond[9], r_cond[49], r_cond[99])
print('維數為10,50,100時的b-Hx_acc矩陣的無窮條件數計算結果：\n', r_B_cond[9], r_B_cond[49], r_B_cond[99])

3.輸出結果

輸入你想計算的Hilbert方陣的最大維數就行，其它交給程式，

請輸入Hilbert方陣的最大維數:100
生成維數從10,20……100的Hilbert矩陣的行范數條件數：
 [35356847610517.12, 6.008376652086652e+18, 8.396589803249062e+18, 9.491653209312077e+19, 1.7763569870536153e+20, 1.9301974218850052e+21, 3.9847310708042826e+19, 1.3450693870678838e+20, 5.444272740462528e+19, 1.3244131088115743e+20]

這里輸出的圖片如下：

這里是第四問的輸出結果：

維數為10,50,100時的x-x_acc計算結果：
 [[-2.54168641e-04]
 [ 2.16242671e-03]
 [-5.54656982e-03]
 [ 5.08880615e-03]
 [ 9.15527344e-04]
 [-4.02832031e-03]
 [ 1.46484375e-03]
 [ 4.88281250e-04]
 [-1.22070312e-04]
 [-6.10351562e-05]] [[ 8.01768149e+02]
 [ 3.33788188e+04]
 [-1.59537467e+06]
 [ 1.98594595e+07]
 [-1.31128704e+08]
·················
 [-4.04700000e+03]
 [ 6.50000000e+01]
 [-2.25000000e+01]
 [ 3.30000000e+01]
 [-7.30000000e+01]] [[ 1.05255071e+04]
 [-1.31934071e+06]
 [ 4.56146227e+07]
 [-7.42843201e+08]
 [ 6.95696228e+09]
 [-4.13027099e+10]
 ················
 [-4.79000000e+02]
 [ 5.12100000e+03]
 [-1.91900000e+03]
 [ 1.77000000e+02]]
 維數為10,50,100時的b-Hx_acc計算結果：
 [[1.27400597e-05]
 [2.15601768e-05]
 [1.73587501e-05]
 [1.43429989e-05]
 [1.23056437e-05]
 [1.08422262e-05]
 [9.72981380e-06]
 [8.84754702e-06]
 [8.12566450e-06]
 [7.52108032e-06]] [[ 13.49201973]
 [  7.51301334]
 [ -4.25849823]
 [-14.49468816]
 [-21.55733783]
 [-26.46828937]
···············
 [-28.3616173 ]
 [-28.05340528]
 [-27.75051982]
 [-27.45291237]] [[-22.02594035]
 [-22.62163018]
 [-20.05848154]
 [-17.70284588]
 [-16.01064382]
 [-14.79343569]
···············
 [ -3.83066178]
 [ -3.79759674]
 [ -3.76500334]
 [ -3.73286444]
 [ -3.70119334]]

維數為10,50,100時的x-x_acc矩陣的無窮條件數計算結果：
 0.00554656982421875 7824513409.0 998313040247.0
維數為10,50,100時的b-Hx_acc矩陣的無窮條件數計算結果：
 2.1560176801216357e-05 38.404672581436365 22.62163018366762

結語

分享出來僅供相互學習，相互探討，如所寫有誤，請多多包涵，希望能相互學習，共同進步，歡迎各位大佬留言評論，