色魔是色魔

几个计算方法上机题目的题解（牛顿迭代/追赶法/三次样条/龙贝格积分）

作者: Semesse
时间: 2019-06-23
分类: Code
评论

炫技就完事了呗，都是 python 写的
***
#### 牛顿迭代法
```python
import math

def sqrtf(c):
    #返回求平方根的f(x)=x^2-c和f'(x)=2x
    return lambda x: x**2-c, lambda x: 2*x

def newton(x0,eps,f,df):
    x=x0-f(x0)/df(x0)
    while(math.fabs(x-x0)>eps):
        x0=x
        x=x0-f(x0)/df(x0)
    return x

n = input("请输入n:")
f,df=sqrtf(n) #求n的平方根
print(newton(1,1e-5,f,df))
```
#### 追赶法解线性方程组
要求系数矩阵是三对角矩阵
```python
def thomas(a, b, c, d):
    n = len(d)
    for i in range(1, n):
        mc = a[i-1]/b[i-1]
        b[i] = b[i] - mc*c[i-1] 
        d[i] = d[i] - mc*d[i-1]
    x = b
    x[-1] = d[-1]/b[-1]
    for i in range(n-2, -1, -1):
        x[i] = (d[i]-c[i]*x[i+1])/b[i]
    return x

if __name__ == '__main__':
    #a:低对角线
    #b:主对角线
    #c:高对角线
    #d:方程右端常数项
    a = list(map(float, input('请输入低对角线元素: ').rstrip().split()))
    b = list(map(float, input('请输入主对角线元素: ').rstrip().split()))
    c = list(map(float, input('请输入高对角线元素: ').rstrip().split()))
    d = list(map(float, input('请输入方程右端常数: ').rstrip().split()))
    print("x = ",thomas(a,b,c,d))
```
#### 三次样条插值
只能算算特定点的值，不能像matlab那样得出系数
```python
from thomas import thomas # 用了上面的追赶法解方程组

#计算二阶差商f[x_i,x_j,x_k]
def f(x,y,i,j,k):
    return y[i]/((x[i]-x[j])*(x[i]-x[k]))+y[j]/((x[j]-x[i])*(x[j]-x[k]))+y[k]/((x[k]-x[i])*(x[k]-x[j]))

def first(x, xs):
    for i in range(len(x)):
        if((xs >= x[i]) and (xs < x[i+1])):
            return i
    else:
        return len(x)

#x,y: list,f(x)=y
#type: 1第一类边界条件 2第二类边界条件
#s0,sn: 边界值
#xs: 要插值的横坐标
def spline(x, y, type, s0, sn, xs):
    n = len(x)-1
    h = [0] + [x[i]-x[i-1] for i in range(1,n+1)]
    # print(n,h)
    if type==1:
        μ = [1] + [h[i]/(h[i]+h[i+1]) for i in range(n)]
        λ = [1-μ[i] for i in range(n)] + [1]
        d = [6*((y[1]-y[0])/(x[1]-x[0])-s0)/h[1]]\
		+ [6*f(x,y,i-1,i,i+1) for i in range(1,n)] + [6*((y[n]-y[n-1])/(x[n]-x[n-1])-sn)/h[n]]
    elif type==2:
        μ = [h[i]/(h[i]+h[i+1]) for i in range(n)]
        λ = [1-μ[i] for i in range(1,n)]
        μ = μ[1:] + [0]
        λ = [0] + λ
        d = [s0] + [6*f(x,y,i-1,i,i+1) for i in range(1,n)] + [sn]
    M = thomas(μ, [2 for _ in range(n+1)], λ, d)
    # print(μ, [2 for _ in range(n+1)], λ, d, M)
    for i in range(len(xs)):
        j = first(x,xs[i])
        print(xs[i],j)
        xs[i] = (x[j+1]-xs[i])**3*M[j]/(h[j+1]*6) + (xs[i]-x[j])**3*M[j+1]/(h[j+1]*6)\
				+ (y[j]-M[j]*h[j+1]*h[j+1]/6)*(x[j+1]-xs[i])/h[j+1] + (y[j+1]-M[j+1]*h[j+1]*h[j+1]/6)*(xs[i]-x[j])/h[j+1]
    return xs
x = list(map(float, input('请输入x坐标(递增)，用空格隔开：').rstrip().split()))
y = list(map(float, input('请输入y坐标，用空格隔开：').rstrip().split()))
type = int(input('请输入边界条件类型(1:夹持边界条件,2:自然边界条件)：'))
s0,sn = list(map(float, input('请输入边界值s0和sn：').rstrip().split()))
xs = list(map(float, input('请输入要插值的x坐标，用空格隔开：').rstrip().split()))
print(spline(x,y,type,s0,sn,xs))
```
#### 龙贝格积分
```python
import numpy as np

#f:待求积函数
#a,b:区间端点
#k:分半次数
def integrate(f, a, b, k):
    t = np.zeros((k, k), dtype=np.float64)
    pow4 = 4 ** np.arange(k, dtype=np.float64) - 1
    h = (b - a)
    t[0, 0] = h * (f(a) + f(b)) / 2
    for j in range(1, k):
        h /= 2
        t[j, 0] = t[j - 1, 0] / 2
        t[j, 0] += h * np.sum(
            f(a + i * h) for i in range(1, 2 ** j + 1, 2)
        )
        for k in range(1, j + 1):
            t[j, k] = t[j, k - 1] + \
                (t[j, k - 1] - t[j - 1, k - 1]) / pow4[k]
    return t[-1, -1]

def f(x):
    return np.exp(x)

print(integrate(f,0,2,5))
```
好了，这就是小编给大家带来的计算方法上机题目的 python 解法。希望大家看完这篇由小编精心整理的内容后，能对相关知识有所了解，用自己的代码写出来。