几个计算方法上机题目的题解（牛顿迭代/追赶法/三次样条/龙贝格积分） 作者: Semesse 时间: 2019-06-23 分类: 千叶 都是 python 写的 <!--more--> *** #### 牛顿迭代法 ```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)) ``` <!--more--> #### 追赶法解线性方程组 要求系数矩阵是三对角矩阵 ```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 解法。希望大家看完这篇由小编精心整理的内容后，能对相关知识有所了解，用自己的代码写出来。 标签: none