|
- # !/usr/bin/env python
- # -*- coding: utf-8 -*-
- """
- @Time : 2020/11/24 10:35
- @Author : Albert Darren
- @Contact : 2563491540@qq.com
- @File : chapter2_3.py
- @Version : Version 1.0.0
- @Description : TODO
- @Created By : PyCharm
- """
- import numpy as np
- import matplotlib.pyplot as plt
- import math
- from matplotlib.font_manager import FontProperties
-
- font_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=12)
-
-
- # 求牛顿n次均差
- def qiujuncha(x, f, n):
- for i in range(1, n):
- for j in range(4, i - 1, -1):
- f[j] = (f[j] - f[j - 1]) / (x[j] - x[j - i])
-
-
- # 根据多项式表达式,x求y
- def function(a, x, n):
- y = 0
- tmp = 1
- for i in range(n):
- y = y + a[i] * tmp
- tmp = tmp * x
- return y
-
-
- # 求x的导数
- def f_dao(x):
- if x == 0:
- return 0
- return 1 / (2 * math.sqrt(x))
-
-
- def draw_pic_niudun(x, s):
- x1 = np.linspace(0, 64, 100)
- plt.plot(x1, lagrange(x1, x, s))
- plt.xlabel(u"x轴", fontproperties=font_set)
- plt.ylabel(u"y轴", fontproperties=font_set)
- plt.show()
-
-
- # 求每一个yk/(xk-x0)(xk-x1)....(xk-xn)
- def qiu_s(s, x, y):
- n = len(x)
- for k in range(n):
- tmp = 1
- for i in range(k):
- tmp = tmp * (x[k] - x[i])
- for i in range(k + 1, n):
- tmp = tmp * (x[k] - x[i])
- s.insert(k, y[k] / tmp)
-
-
- # 根据拉格朗日插值函数求y
- def lagrange(x1, x, s):
- n = len(x)
- result = 0
- for i in range(n):
- tmp = 1
- for j in range(0, i):
- tmp = tmp * (x1 - x[j])
- for j in range(i + 1, n):
- tmp = tmp * (x1 - x[j])
- result += tmp * s[i]
- return result
-
-
- # 第一边界条件下的三次样条插值 求M
- def qiu_m(h, f, o, u, d, x):
- n = len(h)
- o[0] = 1
- u[n] = 1
- d[0] = 6 * (f[0] - f_dao(x[0])) / h[0]
- d[n] = 6 * (f_dao(x[n] - f[n - 1])) / h[n - 1]
- a = []
- for i in range(1, n - 1):
- u[i] = h[i - 1] / (h[i - 1] + h[i])
- for i in range(1, n):
- o[i] = h[i] / (h[i - 1] + h[i])
- for i in range(1, n):
- d[i] = 6 * (f[i + 1] - f[i]) / (h[i - 1] + h[i])
- t = [0 for i in range(n + 1)]
- t[0] = 2
- t[1] = o[0]
- a.append(t)
- for i in range(1, n):
- t = [0 for i in range(n + 1)]
- t[i - 1] = u[i + 1]
- t[i] = 2
- t[i + 1] = o[i + 1]
- a.append(t)
- t = [0 for i in range(n + 1)]
- t[n - 1] = u[n]
- t[n] = 2
- a.append(t)
- tmp = np.linalg.solve(np.array(a), np.array(d))
- m = []
- for i in range(n + 1):
- m.append(tmp[i])
- return m
-
-
- # 根据三次条插值函数求值
- def yangtiao(x1, m, x, y, h, j):
- return m[j] * (x[j + 1] - x1) ** 3 / (6 * h[j]) + m[j + 1] * (x1 - x[j]) ** 3 / (6 * h[j]) + (
- y[j] - m[j] * h[j] ** 2 / 6) * (x[j + 1] - x1) / h[j] + (y[j + 1] - m[j + 1] * h[j] ** 2 / 6) * (
- x1 - x[j]) / h[j]
-
-
- def draw_yangtiao(m, x, y, h, n):
- for i in range(n - 1):
- x1 = np.linspace(x[i], x[i + 1], 10)
- plt.plot(x1, yangtiao(x1, m, x, y, h, i), color='red', label=u"三次样条插值")
- plt.xlabel(u"x轴", fontproperties=font_set)
- plt.ylabel(u"y轴", fontproperties=font_set)
- plt.title(u"三次样条插值", fontproperties=font_set)
- plt.show()
-
-
- def main():
- x = [0, 1, 4, 9, 16, 25, 36, 49, 64]
- y = [0, 1, 2, 3, 4, 5, 6, 7, 8]
- n = len(x)
- s = []
- h = []
- f = y[:]
- for i in range(n - 1):
- h.append(x[i + 1] - x[i])
- u = [0 for n in range(len(h) + 1)]
- d = [0 for n in range(len(h) + 1)]
- o = [0 for n in range(len(h) + 1)]
- qiu_s(s, x, y)
- qiujuncha(x, f, 2)
- m = qiu_m(h, f, o, u, d, x)
- # print(m)
- draw_yangtiao(m, x, y, h, n)
-
-
- main()
|