import numpy as np
import matplotlib.pyplot as plt

x_array = np.arange(1994, 2004, 1)
y_array = np.array([67.052, 68.008, 69.803, 72.024, 73.400, 72.063, 74.669, 74.487, 74.065, 76.777])
m = 10

# 定义插值基函数l
def l(i, x):
    temp0 = 1.0
    temp1 = 1.0
    for k in range(m):
        if k != i:
            temp0 *= (x - x_array[k])
            temp1 *= (x_array[i] - x_array[k])
    return temp0/temp1

# 最终生成的多项式p
def p(x):
    temp = 0.0
    for i in range(m):
        temp += y_array[i]*l(i, x)
    return temp


z1 = np.polyfit(x_array, y_array, 9)
x = np.linspace(1993.5, 2003.5, 100000)
y_val = p(x)
plot1 = plt.plot(x_array, y_array, '*', label='original values')
plot2 = plt.plot(x, y_val, "r")

# 预估2010年的石油量
print(p(2010))
plt.savefig("1-9阶多项式拟合.png")

from scipy.optimize import leastsq
import matplotlib.pyplot as plt
import numpy as np

x_array = np.arange(1994, 2004, 1)
y_array = np.array([67.052, 68.008, 69.803, 72.024, 73.400, 72.063, 74.669, 74.487, 74.065, 76.777])

def aim_func1(p, x):
    """
    一次多项式拟合的方程
    """
    a1, a0 = p
    return a1*x + a0

def aim_func2(p, x):
    """
    二次多项式(曲线)拟合的方程
    """
    a2, a1 , a0 = p
    return a2*(x**2) + a1*x + a0

def error1(p, x, y):
    return aim_func1(p, x)-y

def error2(p, x, y):
    return aim_func2(p, x)-y

p1 = [0, 0]
p2 = [0, 0, 0]
Para1 = leastsq(error1, p1, args=(x_array, y_array))
Para2 = leastsq(error2, p2, args=(x_array, y_array))
a11, a10 = Para1[0]

def f1(x):
    return a11*x + a10
mean1 = 0
mean_square_error1 = 1
for i in range(len(x_array)):
    mean1 += (f1(x_array[i]) - y_array[i])**2
    mean_square_error2 = mean1**0.5
print("拟合为直线时的均方差: " + str(mean_square_error1))
print("估计2010年的值为: " + str(f1(2010)))

a22, a21, a20 = Para2[0]


def f2(x):
    return a22*(x**2) + a21*x + a20
mean2 = 0
mean_square_error2 = 1
for i in range(len(x_array)):
    mean2 += (f2(x_array[i]) - y_array[i])**2
    mean_square_error2 = mean2**0.5
print("拟合曲线为抛物线时的均方差: " + str(mean_square_error2))
print("估计2010年的值为: " + str(f2(2010)))


def three_times_fit():
    # 初始化
    x_array = np.arange(4, 14, 1)
    y_array = np.array([67.052, 68.008, 69.803, 72.024, 73.400, 72.063, 74.669, 74.487, 74.065, 76.777])
    m = len(x_array)
    a = np.zeros((4, 4), dtype='float64')
    b = np.zeros((4, 1))

    def f(i, x):
        """
        基函数span{1,x,x^2,x^3}
        """
        if i == 0:
            return 1
        elif i == 1:
            return x
        elif i == 2:
            return x ** 2
        else:
            return x ** 3

    def compute_eq():
        """
        计算法方程组的两个相关矩阵
        并求出稀疏矩阵res
        :return:
        """
        for i in range(4):
            for j in range(4):
                for x in x_array:
                    a[i][j] += f(i, x) * f(j, x)

        for i in range(4):
            for j in range(len(y_array)):
                b[i] += y_array[j] * f(i, x_array[j])

        a_inv = np.linalg.inv(a)
        res = np.dot(a_inv, b)
        return res

    res = compute_eq()

    def fit_y(x_r):
        x = x_r - 1990
        return res[0] + res[1] * x + res[2] * (x ** 2) + res[3] * (x ** 3)

    x = np.linspace(1994, 2006, 100000)
    plt.plot(x, fit_y(x), color="black", label="y=ax^3+bx^2+cx+d")

    def mean_square_error(x_array, y_array, y):
        mean_error = 0
        for i in range(len(x_array)):
            mean_error += (y_array[i] - y(x_array[i]+1990))**2
            return mean_error**0.5

    mean_square_error3 = mean_square_error(x_array, y_array, fit_y)
    print("拟合三次的均方差是: " + str(float(mean_square_error3)))
    fit2010 = fit_y(2010)
    print("估计2010年的值为: " + str(float(fit2010)))


plt.figure(figsize=(8, 6))
three_times_fit()
plt.scatter(x_array, y_array, color="red", label="Sample Point")
x = np.linspace(1994, 2006, 100000)
y2 = a22*(x*x) + a21*x + a20
plt.plot(x, y2, color="orange", label="y=ax^2+bx+c")
y1 = a11*x + a10
plt.plot(x, y1, color="blue", label="ax+b")

plt.legend()
plt.savefig("2-最小二乘法拟合123次曲线.png")

"""
牛顿迭代法求解方程组x^3+x^2+x-3=0 的根，初始值为x0 = -0.7
迭代7步，并于真值 x*=1 比较
"""

x_true = 1

def f(x):
    """
    y=x^3+x^2+x-3
    :param x:
    :return:y
    """
    return x**3 + x**2 + x - 3


def f_derivatives(x):
    """
    :param x:
    :return: y的导数值
    """
    return 3*(x**2) + 2*x + 1


def f_derivatives2(x):
    """
    :param x:
    :return: y的二阶导数
    """
    return 6*x + 2


def newton_iteration(x):
    """
    牛顿迭代法
    :param x:
    :return:
    """
    return x - f(x) / f_derivatives(x)


if __name__ == '__main__':
    xk = -0.7
    ei_s = []
    for i in range(1, 8):
        print(" "*20 + "第 %s 次迭代" % i + " "*20)
        xk = newton_iteration(xk)
        print("xk的值= " + str(xk))
        ei = abs(xk - x_true)
        print("ei = |xi - x*|= " + str(ei))
        ei_s.append(ei)
        if i != 1:
            print("ei / (last_ei)^2 的值= " + str(ei_s[i-1]/(ei_s[i-2])**2))
        print("")

    print("f(x∗)的二阶导/2(f(x∗)的导数)= " + str(f_derivatives2(x_true)/(2*f_derivatives(x_true))))

import numpy as np
import copy


def gauss(n):
    l = np.zeros((n, n))
    u = np.zeros((n, n))
    a = np.zeros((n, n))
    b = np.ones((n, 1))

    # 初始化系数矩阵A
    for i in range(n):
        for j in range(n):
            a[i][j] = 1/(i+1+j+1-1)

    # 用Doolittle公式求解l，和u的系数
    for j in range(n):
        u[0][j] = a[0][j]
    for j in range(1, n):
        l[j][0] = a[j][0]/u[0][0]
    for d in range(n):
        l[d][d] = 1

    for i in range(1, n):
        for j in range(i, n):
            temp = 0
            for k in range(i):
                temp += l[i][k] * u[k][j]
            u[i][j] = a[i][j] - temp

        for j in range(i+1, n):
            temp = 0
            for k in range(i):
                temp += l[j][k] * u[k][i]
            l[j][i] = (a[j][i]-temp)/u[i][i]

    u_inv = np.linalg.inv(u)
    l_inv = np.linalg.inv(l)
    x = np.dot(u_inv, l_inv)
    x = np.dot(x, b)
    return x


def jacobi(n):
    a = np.zeros((n, n))
    b = np.ones((n, 1))
    # 系数矩阵A初始化
    for i in range(n):
        for j in range(n):
            a[i][j] = 1/(i+1+j+1-1)
    # 初始解x0 = (0,0,0...共n个0)
    x1 = np.array([[0.0], [0.0], [0.0], [0.0]])
    for _ in range(200):
        x = copy.copy(x1)
        for i in range(n):
            temp = 0
            for j in range(n):
                if i != j:
                    temp += a[i][j]*x[j][0]
            x1[i][0] = (b[i][0] - temp)/a[i][i]
    return x1


def gauss_seidel(n):
    a = np.zeros((n, n))
    b = np.ones((n, 1))
    # 系数矩阵A初始化
    for i in range(n):
        for j in range(n):
            a[i][j] = 1/(i+1+j+1-1)
    # 初始解x0 = (0,0,0...共n个0)
    x = np.ones((n, 1))
    for _ in range(7000):
        for i in range(n):
            temp = 0
            for j in range(n):
                temp += a[i][j]*x[j][0]
            x[i][0] = x[i][0] + (b[i][0] - temp)/a[i][i]
    return x


def inner_product(a, b):
    temp = 0
    for i in range(len(a)):
        temp += a[i]*b[i]
    return temp


def c_g(n):
    a = np.zeros((n, n))
    b = np.ones((n, 1))
    # 系数矩阵A初始化
    for i in range(n):
        for j in range(n):
            a[i][j] = 1 / (i + 1 + j + 1 - 1)
    # 初始解x0 = (0,0,0...共n个0)
    r0 = np.zeros((n, 1))
    x = np.ones((n, 1))
    r = b - np.dot(a, x)
    p = r
    time = 0
    while True:
        if (r == r0).all():
            print("C-G法迭代了 " + str(time) + " 次 ")
            return x
        else:
            ap = np.dot(a, p)
            if inner_product(p, ap) == 0:
                print("C-G法迭代了 " + str(time) + " 次 ")
                return x
            else:
                alpha = inner_product(r, r)/(inner_product(p, ap))
                x1 = x + alpha*p
                r1 = r - alpha*ap
                beat = inner_product(r1, r1)/inner_product(r, r)
                p1 = r1 + beat*p
        r = copy.copy(r1)
        x = copy.copy(x1)
        p = copy.copy(p1)
        time += 1


if __name__ == '__main__':
    print(gauss(4))
    print(gauss(8))
    # jacobi(4)
    # jacobi(8)
    print(gauss_seidel(4))
    print(gauss_seidel(8))
    print(c_g(4))
    print(c_g(8))

import numpy as np
from sympy import *


# 初始化矩阵, n为阶数, Am=g即为要求的最终形式, lambda1, miu为中间变量
n = 4
x = np.array([0, 1, 2, 3, 4])
y = np.array([1, 3, 3, 4, 2])
A = np.zeros((5, 5))
lambda1 = np.zeros(5)
miu = np.zeros(5)
h = np.zeros(5)
g = np.zeros(5)

# 根据课本180/181 页公式初始化g, A
for i in range(1, n):
    lambda1[i] = 0.5
    miu[i] = 0.5

for i in range(n):
    h[i] = 1

g[0] = 6
g[4] = -6
for i in range(1, n):
    g[i] = 3*(miu[i]*(y[i+1]-y[i])/h[i] + lambda1[i]*(y[i]-y[i-1])/h[i-1])

A[0][0] = 2
A[0][1] = 1
A[4][3] = 1
A[4][4] = 2
for i in range(1, n):
    A[i][i] = 2
    A[i][i-1] = lambda1[i]
    A[i][i+1] = miu[i]

# 解出m
A_inv = np.linalg.inv(A)
m = np.dot(A_inv, g)

# 确定n个区间内的函数
xx = Symbol('x')
for i in range(n):
    s = ((h[i]+2*(xx - x[i+1]))/(h[i]**3))*((xx - x[i+1])**2)*y[i] + \
        ((h[i]-2*(xx - x[i+1]))/(h[i]**3))*((xx - x[i])**2)*y[i+1] + \
        (xx - x[i])*((xx-x[i+1])**2)*m[i]/(h[i]**2) + (xx - x[i+1])*((xx-x[i])**2)*m[i+1]/(h[i]**2)
    print(s)

"""
用不同的初始值,并使用弦截法,来计算方程组y=x^3 + 2x^2 + 10x - 100的根
"""

import matplotlib.pyplot as plt
import numpy as np


plt.figure(figsize=(8, 6))
x = np.linspace(-2, 4, 1000000)
y = x**3 + 2*(x**2) + 10*x - 100
plt.plot(x, y, color="orange", label="y=x^3 + 2x^2 + 10x - 100")
plt.legend()
plt.savefig("6-弦截法求方程的实根-画图判别个数.png")

"""
经画图判别上述方程只有一个实根(且实根在3.5附近)，可用简单的弦截法计算。
"""


def f(x):
    return x**3 + 2*(x**2) + 10*x - 100


def secant_method(x0, x1):
    return x1 - f(x1)/(f(x1) - f(x0))*(x1 - x0)


def re_secant_method(x0, x1):
    time = 0
    while True:
        temp = secant_method(x0, x1)
        x0 = x1
        x1 = temp
        time += 1
        if abs(f(x1)) - 0 < 0.00001:
            print("迭代了 " + str(time) + " 次")
            print(x1)
            break


if __name__ == '__main__':
    print("当初始值选5, 4时 ")
    re_secant_method(5, 4)
    print("当初始值选100, 99时 ")
    re_secant_method(100, 99)
    print("当初始值选取-100, 50时 ")
    re_secant_method(-100, 50)
    print("当初始值选取100000, 200000时 ")
    re_secant_method(100000, 200000)

import numpy as np


def f(x):
    return (np.exp(3*x))*(np.cos(np.pi*x))


def composite_trapezoidal(n):
    """
    分成n份的复合梯形公式,区间为(0,2pi)
    :param n:
    :return:
    """
    a = 0
    b = 2*np.pi
    h = (b-a)/n
    temp = 0
    for i in range(1, n):
        temp += f(h*i)
    return ((b-a)/(2*n))*(f(a)+f(b)+2*temp)


def composite_simpson(n):
    """
    复合Simpson
    """
    a = 0
    b = np.pi*2
    h = (b-a)/n
    temp1 = 0
    temp2 = 0
    for i in range(1, n):
        temp1 += f(h*i)
    for i in range(0, n):
        temp2 += f((0.5+i)*h)
    return ((b-a)/(6*n))*(f(a)+2*temp1+4*temp2+f(b))


if __name__ == '__main__':
    test_n = [50, 100, 200, 500, 1000]
    for i in range(len(test_n)):
        print("取等分区间n= " + str(test_n[i]) + " 时的积分结果")
        print(composite_trapezoidal(test_n[i]))
        print(composite_simpson(test_n[i]))

"""
本例采用Gauss-legendre方法构造Gauss积分
Gauss求积分的步骤
1. 求积区间和权函数构造n+1次正交多项式
2. 解出正交多项式的 n+1个零点作为插值节点
3. 求解n次代数精度得到的n+1个线性方程
"""
import math

# 被积函数1
def fun1(x):
    return x**2/((1-x**2)**0.5)

# 被积函数2
def fun2(x):
    return math.sin(x)/x


def main(fun, a, b):
    fun = fun
    a = a
    b = b
    # 构造对应的legendre多项式的求积节点
    g_l_2 = {0.5773502692: 1}
    g_l_3 = {0.7745966692: 0.555555556, 0: 0.8888888889}
    g_l_5 = {0.9061798459: 0.2369268851, 0.5384693101: 0.4786286705, 0: 0.5688888889}

    items = ['g_l_2', 'g_l_3', 'g_l_5']

    for i in range(len(items)):
        gauss_sum = 0.0
        functions = [g_l_2, g_l_3, g_l_5]
        for key, value in functions[i].items():
            gauss_sum += fun(((b - a) * key + a + b) / 2) * value
            if key > 0:
                gauss_sum += fun(((a - b) * key + a + b) / 2) * value
        gauss_sum = gauss_sum * (b - a) / 2
        print(items[i], ":   ", gauss_sum)


if __name__ == '__main__':
    print("计算函数1的结果 ")
    main(fun1, -1, 1)
    print("计算函数2的结果 ")
    main(fun2, 0, math.pi*0.5)

"""
y'=y(cost)
y(0)=1
"""
import math
import matplotlib.pyplot as plt
import numpy as np


# 定义y的导数
def y_derivatives(y, t):
    return y*math.cos(t)


# euler法求微分方程
def euler():
    h = [0.1, 0.01, 0.001]
    for h1 in h:
        N = 1/h1
        u = [0 for _ in range(int(N+1))]
        t = [0 for _ in range(int(N+1))]
        u[0] = math.cos(t[0])*1
        for i in range(0, int(N)):
            t[i] = i*h1
            u[i+1] = u[i] + h1*y_derivatives(u[i], t[i])
        plt.scatter(1, u[int(N)], color="red", label="euler")
        print("当步长设置为: " + str(h1) + " Euler法的结果如下")
        print(u[int(N)])


# 改进的euler方法求解微分方程
def euler_plus():
    h = [0.1, 0.01, 0.001]
    for h1 in h:
        N = 1 / h1
        u = [0 for _ in range(int(N + 1))]
        t = [0 for _ in range(int(N + 1))]
        u[0] = 1

        for i in range(0, int(N+1)):
            t[i] = i*h1

        for i in range(0, int(N)):
            u_temp = u[i] + h1*y_derivatives(u[i], t[i])
            u[i+1] = u[i] + (h1*0.5)*(y_derivatives(u[i], t[i])+y_derivatives(u_temp, t[i+1]))

        plt.scatter(1, u[int(N)], color="orange", label="euler_plus")
        print("当步长设置为: " + str(h1) + " 改进的Euler法的结果如下")
        print(u[int(N)])


# 四阶runge_kutta法
def runge_kutta():
    h = [0.1, 0.01, 0.001]
    for h1 in h:
        N = int(1 / h1)
        u = [0 for _ in range(int(N + 1))]
        t = [0 for _ in range(int(N + 1))]
        u[0] = 1
        for i in range(N):
            t[i] = i * h1
            k1 = y_derivatives(u[i], t[i])
            k2 = y_derivatives(u[i]+0.5*h1*k1, t[i]+0.5*h1)
            k3 = y_derivatives(u[i]+0.5*h1*k2, t[i]+0.5*h1)
            k4 = y_derivatives(u[i]+h1*k3, t[i]+h1)
            u[i+1] = u[i]+(h1/6)*(k1+2*k2+2*k3+k4)

        print("当步长设置为: " + str(h1) + " runge_kutta法的结果如下")
        print(u[int(N)])
        plt.scatter(1, u[int(N)], color="blue", label="runge_kutta")


if __name__ == '__main__':
    plt.figure(figsize=(8, 6))
    euler()
    euler_plus()
    runge_kutta()

    print("精确解： ", math.exp(math.sin(1)))
    # 绘图
    x = np.linspace(0, 1, 1000)
    yy = np.exp(np.sin(x))
    plt.plot(x, yy, color="black", label="y=e^sin(t)")
    plt.legend()
    plt.savefig("9-微分方程近似解.png")