跳到主要内容

数值计算期末复习

实验一​

求近似值​

import math


def case_one(num=8):
    I0 = 1 - math.exp(-1.0)
    # I0= input()
    # I0 = float(I0) 读入一个float形式的I0
    # I0 =1-e^-1
    I0 = float("{0:.8f}".format(I0))
    print(I0)
    for k in range(1, num + 1):
        I1 = 1 - k * I0
        I0 = I1
        I0 = float("{0:.8f}".format(I0))
        print("%d, %.8f" % (k, I0))

def case_two(num=8):
    I1 = 0.1124
    for k in range(num, 0, -1):
        I0 = (1 - I1) / k
        I1 = I0
        I1 = float("{0:.4f}".format(I1))
        print("%d, %.4f" % (k - 1, I1))


if __name__ == "__main__":
    chos1 = input()
    chos = int(chos1)
    if chos1 == "1":
        case_one()
    elif chos1 == "2":
        case_two()

第一题读入一个数值 choose 然后选择，第一个函数的推导公式是Tn+1=1−nTnT_{n+1}=1-nT_n，第二个函数的推导公式是Tn=1−Tn−1nT_n=\frac{1-T_{n-1}}{n}，然后输出结果其中中间的过程保留小数通过 format 函数保留小数点后 8 位，第二个函数保留小数点后 4 位。

泰勒展开求解 sin(x)​

import math


def solve_sin(x, change=1e-3):
    i = 1
    # i 是第i项
    sum = 0
    term = x
    while abs(term) > change:
        sum += term
        term *= -x * x / (2 * i) / (2 * i + 1)
        # term = -x * x / (2 * i) / (2 * i + 1) * term
        i += 1
    sum += term
    return round(sum, 8)


def solve_cos(x, change=1e-3):
    i = 1
    # i 是第i项
    sum = 0
    term = 1
    while abs(term) > change:
        sum += term
        term *= -x * x / (2 * i) / (2 * i - 1)
        i += 1
    sum += term
    return round(sum, 8)


x = float(input())
x = math.fmod(x, 2 * math.pi)

print(format(solve_sin(x), ".8f"))
print(format(solve_cos(x), ".8f"))

题目是泰勒展开求解 sin(x)和 cos(x)，然后读入一个数值 x，然后通过 format 函数保留小数点后 8 位，然后输出结果。 format 函数的用法是 format(数值，".nf")，其中 n 是保留小数点后 n 位。

实验二​

二分法求解方程​

# 二分法
import numpy as np
import math



def func(x):
    return np.exp(x) + 10 * x - 2  # 定义函数e


def dichotomy(a, b, tol):
    mid = (a + b) / 2  # 初始化中点
    y1 = func(a)
    y2 = func(b)
    if y1 * y2 > 0:
        return None  # 输出失败信息
    k = 0
    while b - a >= tol:  # b - a < e
        mid = (a + b) / 2
        y = func(mid)
        # print(k, a, b, func(a), func(b), func(mid), mid)
        if abs(y) < 1e-8:  # y=0, 满足精度要求
            break
        if y1 * y < 0:
            b = mid
        else:
            a = mid
            y1 = y
        k = k + 1

    mid = (a + b) / 2
    print("%.5f, %.9f" % (mid, func(mid)))

    return mid


if __name__ == "__main__":
    ## sample input: -1 1 0.5e-3
    a, b, tol = map(float, input().split())
    dichotomy(a, b, tol)

标准的二分求解注意要在最后一次二分之后加上去，其中 total 是精度 func 是定义函数的地方，然后读入一个数值 a,b,tol，然后输出结果。

不动点法求迭代方程​

# 不动点迭代
import math


def ff(x):
    return math.exp(x) + 10 * x - 2


def ff0(x):
    return (2 - math.exp(x)) / 10


def f0(x):
    return x**3 - x - 1.0


def gf0(x, case=1):
    if case == 0:
        return (x + 1.0) ** (1 / 3.0)
    else:
        return x**3 - 1.0


def f1(x):
    return x - x ** (1 / 3.0) - 2


def gf1(x, case=2):
    if case == 0:
        return x ** (1 / 3.0) + 2
    elif case == 1:
        return (x - 2.0) ** 3
    else:
        return (6 + 2 * x ** (1 / 3.0)) / (3 - x ** (-2.0 / 3))


def f2(x):
    return x - math.exp(-x)


def gf2(x):
    return math.exp(-x)


def fixed_point(x0, tol, func, gfunc, N=100):
    for k in range(N):
        # print(k, x0, func(x0))
        # 输出的是迭代次数，当前迭代值，当前迭代值的函数值
        x1 = gfunc(x0)
        if abs(x1 - x0) < tol:
            x0 = x1
            print("%.5f, %.9f" % (x1, func(x1)))
            return 0
        x0 = x1
    print("达到最大迭代次数！")
    return 1


if __name__ == "__main__":
    # x0 = 1.5
    # tol = 1e-5
    # fixed_point(x0, tol, func=f0, gfunc=gf0)

    # x0 = 3.5
    # tol = 1e-8
    # fixed_point(x0, tol, func=f1, gfunc=gf1)

    # x0 = 0.5
    # tol = 1e-5
    # fixed_point(x0, tol, func=f2, gfunc=gf2)
    x0 = 0
    total = 0.5e-3
    fixed_point(x0, total, func=ff, gfunc=ff0)

牛顿迭代法​

# 牛顿迭代法
import numpy as np

# import matplotlib.pyplot as plt


def fun0(x):
    return x**3 - x - 1
    # return x - np.cos(x)


def gfun0(x):
    return 3 * x**2 - 1
    # return 1 + np.sin(x)


def ff(x):
    return np.exp(x) + 10 * x - 2


def ff1(x):
    return np.exp(x) + 10


# def draw_line(x0, x1, func, num=100):
#     xs = np.linspace(x0, x1, num)
#     ys = func(xs)
#     plt.figure()
#     plt.plot(xs, ys, "r")
#     plt.show()


def newton_iteration(x0, tol, gfunc, func, N=100):
    for k in range(N):
        # print("k=%d, x=%.6f, f(x)=%.6f" % (k, x0, func(x0)))
        f0, g0 = func(x0), gfunc(x0)
        if abs(g0) < 1e-9:
            print("梯度为零！")
            return -1
        x1 = x0 - f0 / g0
        if abs(x1 - x0) < tol:
            # print("k=%d, x=%.6f, f(x)=%.6f" % (k + 1, x1, func(x1)))
            print("%.8f" % (x1))
            return 0
        x0 = x1
    print("达到最大迭代次数")
    return 1


if __name__ == "__main__":
    # x0, tol = 1, 1e-5  # x0 = 1, 0, pi/2
    x0, tol = map(float, input().split())
    newton_iteration(x0, tol, func=ff, gfunc=ff1)
    # draw_line(0, np.pi / 2, func=fun0)

这里的 ff1 是由 ff 求导获得的需要注意，然后达到精度范围就输出了，然后读入一个数值 x0,tol，然后输出结果。

牛顿下山法​

import numpy as np


def ff(x):
    return np.exp(x) + 10 * x - 2


def ff1(x):
    return np.exp(x) + 10


def newton_downhill(x0, tol, gfunc, func, N=100):
    for k in range(N):
        # print("k=%d, x=%.6f, f(x)=%.6f" % (k, x0, func(x0)))
        f0, g0 = func(x0), gfunc(x0)
        if abs(g0) < 1e-9:
            print("梯度为零！")
            return -1
        x1 = x0 - f0 / g0
        # 线性搜索
        lam = 1
        for i in range(N):
            x1_ = lam * x1 + (1 - lam) * x0
            #这一步是为了寻找下降的点
            if abs(func(x1_)) < abs(func(x0)):  # 下降条件
                x1 = x1_
                break
            print("%d, %.6f, %.6f, %.6f" % (k, lam, x1_, func(x1_)))
            lam *= 0.5
            pass
        if abs(x1 - x0) < tol:
            print("%.8f, %.6f" % (x1, func(x1)))
            return 0
        x0 = x1
        pass  # pass 是空语句，是为了保持程序结构的完整性
    print("达到最大迭代次数")
    return 1


x0, tol = map(float, input().split())
newton_downhill(x0, tol, func=ff, gfunc=ff1)

牛顿下山法和牛顿迭代法的区别在于牛顿下山法是在牛顿迭代法的基础上加上了线性搜索，然后读入一个数值 x0,tol，然后输出结果。

弦截法​

import numpy as np


def ff(x):
    return np.exp(x) + 10 * x - 2


def ff1(x):
    return np.exp(x) + 10


def secant(x0, tol, gfunc, func, N=100, digit=6):
    fmt = "k=%d, x=%." + str(digit) + "f, f(x)=%." + str(digit) + "f"
    f0, g0 = func(x0), gfunc(x0)
    # print(fmt % (0, x0, f0))
    x1 = x0 - f0 / g0
    f1 = func(x1)
    for k in range(1, N):
        # print(fmt % (k, x1, func(x1)))
        if abs(f1 - f0) < 1e-9:
            print("梯度为零！")
            return -1
        x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
        if abs(x2 - x1) < tol:
            print("%.8f, %.6f" % (x2, func(x2)))
            return 0
        x0, f0 = x1, f1
        x1, f1 = x2, func(x2)
    print("达到最大迭代次数")
    return 1


xo, tol = map(float, input().split())
secant(xo, tol, func=ff, gfunc=ff1)

实验三​

艾特肯法求方程的根​

import numpy as np


# x=ax^3+b
def fun(a: float, b: float, x: float):
    return a * x**3 + b


def aitken(a: float, b: float, x: float):
    lastx = x
    y = x + 100
    z = x + 100
    while True:
        lastx = x
        y = fun(a, b, x)
        z = fun(a, b, y)
        x = x - (y - x) ** 2 / (z - 2 * y + x)  # aitken 的迭代公式
        # print(f"debug : {y:.6f}, {z:.6f}, {x:.6f}")
        if abs(lastx - x) < 1e-5:
            break
            #这里有个精度问题
    print(f"{y:.6f}, {z:.6f}, {x:.6f}")


a, b, x0 = map(float, input().split())
aitken(a, b, x0)

输入的是 a,b,x0，然后输出结果。

求解实部虚部​

def newton(f, df, x0, tol=1e-5, num=1000):
    x = x0
    for i in range(num):
        if df(x) == 0:
            return "error"
        x_new = x - f(x) / df(x)
        if abs(x_new.real - x.real) < tol and abs(x_new.imag - x.imag) < tol:
            x_new = x - f(x) / df(x)
            return "{:.8f} {:.8f}i".format(x_new.real, x_new.imag)
        x = x_new
    return "error"


f = lambda x: x**4 + 2 * x**3 - x - 1
df = lambda x: 4 * x**3 + 6 * x**2 - 1

x0, xi = map(float, input().split())
x0 = complex(x0, xi)
print(newton(f, df, x0))

输入的是实部虚部，然后输出结果。

列主元消元法求解线性方程组​

# 列主元消元法求解线性方程组
import numpy as np

def print_matrix(a):
    for i in a:
        for j in i:
            print(f"{j:.5f}", end=" ")
        print()
    print()

# column pivot gauss
def gauss(a, b):
    n = b.size

    for i in range(n):
        k = i
        for j in range(i + 1, n):
            if abs(a[j, i]) > abs(a[k, i]):
                k = j
        if k != i:
            a[[i, k]] = a[[k, i]]
            b[[i, k]] = b[[k, i]]
        # print_matrix(a)

        for j in range(i + 1, n):
            c = a[j, i] / a[i, i]
            a[j] -= c * a[i]
            b[j] -= c * b[i]

    # print(a)
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - np.dot(a[i, i + 1:], x[i + 1:])) / a[i, i]

    return x


if __name__ == '__main__':
    s = list(map(float, input().split()))
    a = np.array(s)
    s = list(map(float, input().split()))
    b = np.array(s)
    a = a.reshape(b.size, -1)

    x = gauss(a, b)

    for i in x:
        print(f"{i:.5f}")

这里求解的是 Ax=b 的解，输入的是矩阵 A 和 b，然后输出结果。

求矩阵的逆​

import numpy as np
def solve(A):
    A = np.array(A, dtype=float)
    n = A.shape[0]
    attend = np.hstack([A, np.eye(n)])
    for i in range(n):
        max_row = np.argmax(np.abs(attend[i:, i])) + i
        attend[[i, max_row]] = attend[[max_row, i]]
        attend[i] = attend[i] / attend[i, i]
        for j in range(n):
            if i != j:
                attend[j] -= attend[i] * attend[j, i]
    A_inv = attend[:, n:]
    return np.round(A_inv, 5)

def solve2():
    A=input().strip().split()
    size=int(len(A)**0.5)
    A=[list(map(float,A[i*size:(i+1)*size])) for i in range(size)]
    A_1=np.linalg.inv(A)
    for(row) in A_1:
        for x in row:
            print(f"{x:.5f}",end=" ")
        print()

def main():
    A = input().strip().split()
    size = int(len(A) ** 0.5)
    A = [list(map(float, A[i*size:(i+1)*size])) for i in range(size)]
    A_1 = solve(A)
    for row in A_1:
        for x in row:
            print(f"{x:.5f}", end=" ")
        print()



if __name__ == "__main__":
    main()

LU 分解法求解线性方程组​

# LU分解法求解线性方程组
# LU分解法求解线性方程组
import math
import numpy as np


def LU(A):
    n = A.shape[0]
    L = np.identity(n)
    U = np.zeros((n, n))
    for i in range(n):
        L[i, i] = 1
        for j in range(i, n):
            U[i, j] = A[i, j] - np.dot(L[i, :i], U[:i, j])
        for j in range(i + 1, n):
            L[j, i] = (A[j, i] - np.dot(L[j, :i], U[:i, i])) / U[i, i]
    return L, U


def solve(L, U, b):
    n = b.shape[0]
    y = np.zeros(n)
    x = np.zeros(n)
    # 解下三角矩阵方程 Ly = b
    for i in range(n):
        y[i] = b[i] - np.dot(L[i, :i], y[:i])
    # 解上三角矩阵方程 Ux = y
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - np.dot(U[i, i + 1 :], x[i + 1 :])) / U[i, i]
    return x, y


def print_mat(A):
    n = A.shape[0]
    for i in range(n):
        for j in range(n):
            print(f"{A[i][j]:.5f}", end=" ")
        print()


def print_vec(b):
    n = b.shape[0]
    for i in range(n):
        print(f"{b[i]:.5f}")


if __name__ == "__main__":
    # 输入矩阵
    A = list(map(float, input().split()))
    A = np.array(A)

    b = list(map(float, input().split()))
    b = np.array(b)

    n = int(math.sqrt(A.shape[0]))
    A = A.reshape(n, n)

    # LU 分解
    L, U = LU(A)
    # 求解
    x, y = solve(L, U, b)
    print("L:")
    print_mat(L)
    print("U:")
    print_mat(U)
    print("y:")
    print_vec(y)
    print("x:")
    print_vec(x)

追赶法求解三对角线性方程组​

# 追赶法求解线性方程组
import numpy as np

# 三对角矩阵追赶法LU分解
def LU(A):
    # 获取矩阵大小
    n = A.shape[0]

    # 初始化 L 和 U 矩阵
    L = np.identity(n)
    U = np.zeros((n, n))

    # 进行分解
    for i in range(n):
        if i == 0:
            # 第一行特殊处理
            U[i][i] = A[i][i]
            U[i][i+1] = A[i][i+1]
        elif i == n-1:
            # 最后一行特殊处理
            L[i][i-1] = A[i][i-1] / U[i-1][i-1]
            U[i][i] = A[i][i] - L[i][i-1] * A[i-1][i]
        else:
            # 中间行通用处理
            L[i][i-1] = A[i][i-1] / U[i-1][i-1]
            U[i][i] = A[i][i] - L[i][i-1] * A[i-1][i]
            U[i][i+1] = A[i][i+1]

    return L, U

def solve(L, U, b):
    # 解下三角矩阵方程 Ly = b
    n = len(L)
    y = np.zeros(n)
    y[0] = b[0]
    for i in range(1, n):
        y[i] = b[i] - L[i][i-1] * y[i-1]

    # 解上三角矩阵方程 Ux = y
    x = np.zeros(n)
    x[n-1] = y[n-1] / U[n-1][n-1]
    for i in range(n-2, -1, -1):
        x[i] = (y[i] - U[i][i+1] * x[i+1]) / U[i][i]

    return x, y

def print_mat(A):
    n = A.shape[0]
    for i in range(n):
        for j in range(n):
            print(f"{A[i][j]:.5f}", end = ' ')
        print()

def print_vec(b):
    n = b.shape[0]
    for i in range(n):
        print(f"{b[i]:.5f}")

if __name__ == '__main__':
    # 输入矩阵
    A = list(map(float, input().split()))
    A = np.array(A)
    b = list(map(float, input().split()))
    b = np.array(b)
    n = b.shape[0]
    A = A.reshape(n, n)
    # LU 分解
    L, U = LU(A)
    # 求解
    x, y = solve(L, U, b)
    print("L:")
    print_mat(L)
    print("U:")
    print_mat(U)
    print("y:")
    print_vec(y)
    print("x:")
    print_vec(x)

平方根法求解线性方程组​

# 平方根法求解线性方程组
import numpy as np

def cholesky(A):
    n = A.shape[0]
    L = np.zeros((n, n))
    for j in range(n):
        s = sum(L[j, k]*L[j, k] for k in range(j))
        L[j, j] = np.sqrt(A[j, j] - s)

        for i in range(j, n):
            t = sum(L[i, k]*L[j, k] for k in range(j))
            L[i, j] = (A[i, j] - t) / L[j, j]

    return L

def solve(L, b):
    n = b.shape[0]
    y = np.zeros(n)
    x = np.zeros(n)

    y[0] = b[0] / L[0][0]
    for i in range(1, n):
        y[i] = (b[i] - sum(L[i][k]*y[k] for k in range(i))) / L[i][i]

    x[n-1] = y[n-1] / L[n-1][n-1]
    for i in range(n-2, -1, -1):
        x[i] = (y[i] - sum(L[k][i]*x[k] for k in range(i+1, n))) / L[i][i]

    return x, y

def print_mat(A):
    n = A.shape[0]
    for i in range(n):
        for j in range(n):
            print(f"{A[i][j]:.5f}", end = ' ')
        print()

def print_vec(b):
    n = b.shape[0]
    for i in range(n):
        print(f"{b[i]:.5f}")

if __name__ == '__main__':
    # 输入矩阵
    A = list(map(float, input().split()))
    A = np.array(A)
    b = list(map(float, input().split()))
    b = np.array(b)
    n = b.shape[0]
    A = A.reshape(n, n)

    L = cholesky(A)
    # 求解
    x, y = solve(L, b)

    print("L:")
    print_mat(L)
    print("y:")
    print_vec(y)
    print("x:")
    print_vec(x)

Jacobi 迭代法求解线性方程组​

# Jacobi法求解线性方程组
import numpy as np

def inf_norm(x):
    return max(abs(x))

def jacobi(A, b, x0):
    n = len(b)
    x = x0
    x_new = np.zeros(n, dtype=np.float32)
    for i in range(100):
        for j in range(n):
            sum = 0
            for k in range(n):
                if k != j:
                    sum += A[j, k] * x[k]
            x_new[j] = (b[j] - sum) / A[j, j]
        if inf_norm(x - x_new) < 1e-5:
            x = x_new.copy()
            break
        x = x_new.copy()
    return x

def print_vec(b):
    n = b.shape[0]
    for i in range(n):
        print(f"{b[i]:.5f}")

if __name__ == '__main__':
    # 输入矩阵
    A = list(map(float, input().split()))
    A = np.array(A)
    b = list(map(float, input().split()))
    b = np.array(b)
    n = b.shape[0]
    A = A.reshape(n, n)
    x0 = list(map(float, input().split()))
    x0 = np.array(x0)

    x = jacobi(A, b, x0)
    print("x:")
    print_vec(x)

Gauss-Seidel 迭代法求解线性方程组​

# G-S法求解线性方程组
import numpy as np

def inf_norm(x):
    return max(abs(x))

def GS(A, b, x0):
    n = len(b)
    x = x0.copy()
    for i in range(100):
        for j in range(n):
            sum = 0
            for k in range(n):
                if k != j:
                    sum += A[j, k] * x[k]
            x[j] = (b[j] - sum) / A[j, j]
        if inf_norm(x - x0) < 1e-5:
            break
        x0 = x.copy()
    return x

def print_vec(b):
    n = b.shape[0]
    for i in range(n):
        print(f"{b[i]:.5f}")

if __name__ == '__main__':
    # 输入矩阵
    A = list(map(float, input().split()))
    A = np.array(A)
    b = list(map(float, input().split()))
    b = np.array(b)
    n = b.shape[0]
    A = A.reshape(n, n)
    x0 = list(map(float, input().split()))
    x0 = np.array(x0)

    x = GS(A, b, x0)
    print("x:")
    print_vec(x)

实验四​

拉格朗日插值法​

# 拉格朗日插值
def lagrange(x, y, x0):
    n = len(x)
    s = 0
    for i in range(n):
        t = 1
        for j in range(n):
            if i != j:
                t *= (x0 - x[j]) / (x[i] - x[j])
        s += t * y[i]
    return s

if __name__ == '__main__':
    x = [0.0,0.1,0.195,0.3,0.401,0.5]
    y = [0.39894,0.39695,0.39142,0.38138,0.36812,0.35206]
    x0 = float(input())
    y0 = lagrange(x, y, x0)
    print(f"{y0:.5f}")

牛顿插值法​

# 牛顿型插值多项式
import math

# 差商
def diff_quotient(x, y):
    n = len(x)
    d = [[0] * n for _ in range(n)]
    for i in range(n):
        d[i][0] = y[i]
    for j in range(1, n):
        for i in range(n - j):
            d[i][j] = (d[i + 1][j - 1] - d[i][j - 1]) / (x[i + j] - x[i])
    return d


def calc_diff_quotient(d, x0):
    n = len(d)
    s = d[0][0]
    for i in range(1, n):
        t = 1
        for j in range(i):
            t *= (x0 - x[j])
        s += d[0][i] * t
    return s


# 差分
def diff(x, y):
    n = len(x)
    d = [[0] * n for _ in range(n)]
    for i in range(n):
        d[i][0] = y[i]
    for j in range(1, n):
        for i in range(n - j):
            d[i][j] = d[i + 1][j - 1] - d[i][j - 1]
    return d


def calc_diff(x, d, x0):
    n = len(d)
    s = d[0][0]
    t = 0
    for i in range(0, n):
        if x0 > x[i]:
            t = (x0 - x[i]) / (x[1] - x[0])
            break

    for i in range(1, n):
        tt = 1
        for j in range(i):
            tt *= (t - j)
        s += d[0][i] * tt / math.factorial(i)
    return s


def print_mat(mat):
    for row in mat:
        for e in row:
            print(f"{e:.5f}", end=' ')
        print()


if __name__ == '__main__':
    x = list(map(float, input().split()))
    y = list(map(float, input().split()))
    x0 = float(input())

    d_fq = diff_quotient(x, y)
    y0_fq = calc_diff_quotient(d_fq, x0)
    d_ff = diff(x, y)
    y0_ff = calc_diff(x, d_ff, x0)

    print_mat(d_fq)
    print(f"{y0_fq:.5f}")
    print_mat(d_ff)
    print(f"{y0_ff:.5f}")

三次样条插值​

# 三次样条插值
import numpy as np

# 三对角矩阵追赶法LU分解
def LU(A):
    # 获取矩阵大小
    n = A.shape[0]

    # 初始化 L 和 U 矩阵
    L = np.identity(n)
    U = np.zeros((n, n))

    # 进行分解
    for i in range(n):
        if i == 0:
            # 第一行特殊处理
            U[i][i] = A[i][i]
            U[i][i+1] = A[i][i+1]
        elif i == n-1:
            # 最后一行特殊处理
            L[i][i-1] = A[i][i-1] / U[i-1][i-1]
            U[i][i] = A[i][i] - L[i][i-1] * A[i-1][i]
        else:
            # 中间行通用处理
            L[i][i-1] = A[i][i-1] / U[i-1][i-1]
            U[i][i] = A[i][i] - L[i][i-1] * A[i-1][i]
            U[i][i+1] = A[i][i+1]

    return L, U

def solve(L, U, b):
    # 解下三角矩阵方程 Ly = b
    n = len(L)
    y = np.zeros(n)
    y[0] = b[0]
    for i in range(1, n):
        y[i] = b[i] - L[i][i-1] * y[i-1]

    # 解上三角矩阵方程 Ux = y
    x = np.zeros(n)
    x[n-1] = y[n-1] / U[n-1][n-1]
    for i in range(n-2, -1, -1):
        x[i] = (y[i] - U[i][i+1] * x[i+1]) / U[i][i]

    return x, y

def spline3(x, y, l, r, x0):
    n = len(x)
    h = np.zeros(n)
    for i in range(1, n):
        h[i] = x[i] - x[i - 1]

    mul = np.zeros(n)
    lam = np.zeros(n)
    for i in range(1, n - 1):
        mul[i] = h[i] / (h[i] + h[i + 1])
        lam[i] = 1 - mul[i]

    g = np.zeros(n)
    g[0] = 6 / h[1] * ((y[1] - y[0]) / h[1] - l)
    g[-1] = 6 / h[-1] * (r - (y[-1] - y[-2]) / h[-1])
    for i in range(1, n - 1):
        g[i] = 6 / (h[i] + h[i + 1]) * ((y[i + 1] - y[i]) / h[i + 1] - (y[i] - y[i - 1]) / h[i])

    d = np.zeros(n)
    d[0] = 6 * ((y[1] - y[0]) / h[1] - l) / h[1]
    d[n - 1] = 6 * (r - (y[n - 1] - y[n - 2]) / h[n - 1]) / h[n - 1]

    A = np.identity(n) * 2
    A[0][1] = 1
    A[-1][-2] = 1
    for i in range(1, n - 1):
        A[i][i - 1] = mul[i]
        A[i][i + 1] = lam[i]

    # M = np.linalg.solve(A, g)
    L, U = LU(A)
    M, _ = solve(L, U, g)

    # print(A)
    # print(g)
    # print(M)

    y0 = 0

    for i in range(n - 1):
        if x0 >= x[i] and x0 <= x[i + 1]:
            y0 += M[i] * (x[i + 1] - x0)**3 / (6 * h[i + 1])
            y0 += M[i + 1] * (x0 - x[i])**3 / (6 * h[i + 1])
            y0 += (y[i] - M[i] * h[i + 1]**2 / 6) * (x[i + 1] - x0) / h[i + 1]
            y0 += (y[i + 1] - M[i + 1] * h[i + 1]**2 / 6) * (x0 - x[i]) / h[i + 1]
            break

    return M, y0


def print_vec(b):
    for i in b:
        print(f"{i:.5f}")

if __name__ == '__main__':
    x = list(map(float, input().split()))
    y = list(map(float, input().split()))
    l, r = map(float, input().split())
    x0 = float(input())

    M, y0 = spline3(x, y, l, r, x0)
    print_vec(M)
    print(f"{y0:.5f}")
0%