[计算方法]上机

我思故我在Fighting

已于 2023-02-02 16:01:59 修改

阅读量220

点赞数

文章标签： python numpy

于 2022-11-10 22:57:46 首次发布

本文链接：https://blog.csdn.net/qq_61897765/article/details/127798084

版权

p11上机一

1.1

def f1(x):
    y=x*((x+1)**0.5-x**0.5)
    return y
def f2(x):
    y=x/((x+1)**0.5+x**0.5)
    return y
f_1=0.414213562
f_2=50000

print(f'公式一f(1)近似值{f1(1)},公式一f(10**10)近似值{f1(10**10)}')
print(f'公式二f(1)近似值{f1(1)},公式二f(10**10)近似值{f2(10**10)}')
print(f'公式一f(1)绝对误差{abs(f1(1)-f_1)},公式一f(10**10)绝对误差{abs(f1(10**10)-f_2)}')
print(f'公式二f(1)绝对误差{abs(f2(1)-f_1)},公式二f(10**10)绝对误差{abs(f2(10**10)-f_2)}')
print(f'公式一f(1)相对误差{abs(f1(1)-f_1)/f_1},公式一f(10**10)相对误差{abs(f1(10**10)-f_2)/f_2}')
print(f'公式二f(1)相对误差{abs(f2(1)-f_1)/f_1},公式二f(10**10)相对误差{abs(f2(10**10)-f_2)/f_2}')

1.2

import math

i0=math.log(1.2)
a=i0
for i in range(1,21):
    b=-5*a+1/i
    print(f'I{i}={b}')
    a=b

p29上机二

2.1二分法

def f(x):
    y = x**3 + 4 * x**2 - 10
    return y

def bisection(a, b, c):
    while abs(b - a)/2 >= c:
        s = (a + b) / 2
        if f(a) * f(s) > 0:
            a = s
        else:
            b = s
        print(s)
    s = (a + b) / 2
    print(s)
a, b, c = 1, 2, 0.00025
bisection(a, b, c)

2.2迭代法

import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
def g1(x):
    xk = x**2 - 1
    return xk
def g2(x):
    xk = 1 + 1 / x
    return xk
def g3(x):
    xk = (x + 1)**0.5
    return xk
chu = 1.3
plt.title('初值为0.9下三个迭代函数的比较')
plt.xlabel('迭代次数')
plt.ylabel('不同迭代函数')
wucha = 0.000001
x0=chu
i = 0
y_list1 = [x0]
x_list1= [i]
while 1:
    i=i+1
    res = g1(x0)
    diff = abs(res - x0)
    x0 = res
    y_list1.append(x0)
    x_list1.append(i)
    # print(x0)
    if diff < wucha:
        break
    if i > 10:
        break
x_list1 = np.array(x_list1)
y_list1 = np.array(y_list1)
x0=chu
i = 0
y_list2 = [x0]
x_list2 = [i]
while 1:
    i=i+1
    res = g2(x0)
    diff = abs(res - x0)
    x0 = res
    y_list2.append(x0)
    x_list2.append(i)
    # print(x0)
    if diff < wucha:
        break
    if i > 10:
        break
x_list2 = np.array(x_list2)
y_list2 = np.array(y_list2)
x0=chu
i = 0
y_list3 = [x0]
x_list3 = [i]
while 1:
    i=i+1
    res = g3(x0)
    diff = abs(res - x0)
    x0 = res
    y_list3.append(x0)
    x_list3.append(i)
    # print(x0)
    if diff < wucha:
        break
    if i > 10:
        break
x_list3 = np.array(x_list3)
y_list3 = np.array(y_list3)
plt.plot(x_list1, y_list1,x_list2, y_list2,x_list3, y_list3)
plt.legend(["迭代函数1", "迭代函数2","迭代函数3"], loc='lower left')
plt.show()

(2)从0到3每隔0.5选取一个数作为初值对迭代函数3进行分析

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False


def g(x):
    xk = (x + 1) ** 0.5
    return xk


x0_list = list(np.arange(0, 3, 0.5))
plt.title('迭代函数3')
plt.xlabel('迭代次数')
plt.ylabel('不同初值')
for x0 in x0_list:
    i = 0
    la = x0
    y_list = [x0]
    x_list = [i]
    wucha = 0.000001
    while 1:
        i = i + 1
        res = g(x0)
        diff = abs(res - x0)
        x0 = res
        y_list.append(x0)
        x_list.append(i)
        print(x0)
        if diff < wucha:
            break
    x_list = np.array(x_list)
    y_list = np.array(y_list)
    plt.plot(x_list, y_list, label=f'初值{la}')
    plt.legend()
plt.show()

2.3牛顿法根号a

import matplotlib.pyplot as plt
import numpy as np
def f(x, a):
    y = x ** 2 - a
    return y


def f1(x):
    y = 2 * x
    return y


def diedai(x0, ep, a):
    x1 = x0 - f(x0, a) / f1(x0)
    k = 1
    y_list = [x0]
    x_list = [k]
    while (abs(x1 - x0) >= ep):
        k = k + 1
        x1 = x0
        G = f1(x1)
        x0 = x1 - f(x1, a) / G
        y_list.append(x0)
        x_list.append(k)
        print(x0)
    x_list = np.array(x_list)
    y_list = np.array(y_list)
    plt.plot(x_list, y_list, label=f'x0={y_list[0]}')
    plt.legend()


x0 = 1.5
ep = 0.00000001
a = 2
diedai(x0, ep, a)

p69 上机三 3.3 随机数据的生成

A = np.zeros((20,20))
for i in range(20):
    for j in range(20):
        if abs(j-i)==2:
            A[i][j]=-0.25
        elif abs(j-i)==1:
            A[i][j]=-0.5
        elif i==j:
            A[i][j]=3
B=np.random.randn(20)*10
x0=np.random.randn(20)*1.5
x=np.zeros((20))
esp=10e-5

雅可比迭代法：

def jacobi(A, B, x0, x, eps):  # A为系数矩阵 B为右端向量 x0为初始向量 x为迭代向量
    times = 0

    while True:
        for i in range(3):
            temp = 0
            for j in range(3):
                if i != j:
                    temp += x0[j] * A[i][j]
            x[i] = (B[i] - temp) / A[i][i]
        calTemp = max(abs(x - x0))
        times += 1
        if calTemp < eps:
            break
        else:
            x0 = x.copy()
  print(times)
    print(x)
    # ylist1.append(times)
    cishu.append(times)
if __name__ == '__main__':
    jacobi(A, B, x0, x,eps)

高斯-赛德尔：

import numpy as np
import matplotlib.pyplot as plt
A = np.zeros((20,20))
for i in range(20):
    for j in range(20):
        if abs(j-i)==2:
            A[i][j]=-0.25
        elif abs(j-i)==1:
            A[i][j]=-0.5
        elif i==j:
            A[i][j]=3
B=np.random.randn(20)*10
x0=np.random.randn(20)*1.5
x=np.zeros((20))
eps=10e-5
#y1data=[]
#y2data=[]
C=A
def jacobi(A, B, x0, x, eps): # A为系数矩阵 B为右端向量 x0为初始向量 x为迭代向量
    times = 0
    while True:
        for i in range(3):
            temp = 0
            for j in range(3):
                if i != j:
                    temp += x0[j] * A[i][j]
            x[i] = (B[i] - temp) / A[i][i]
        calTemp = max(abs(x - x0))
        times += 1
        if calTemp < eps:
            break
        else:
            x0 = x.copy()
    print(times)
    print(x)
    y1data.append(times)

def gass(A, B, x0, x,eps):
    times = 0
    while True:
        for i in range(3):
            temp = 0
            tempx = x0.copy()
            for j in range(3):
                if i != j:
                    temp += x0[j] * A[i][j]
            x[i] = (B[i] - temp) / A[i][i]
            x0[i] = x[i].copy()
        calTemp = max(abs(x - tempx))
        times += 1
        if calTemp < eps:
            break
        else:
            x0 = x.copy()
    print(times)
    print(x)
    y2data.append(times)
if __name__ == '__main__':
    #第一问
    y1data = []
    y2data = []

    for i in range(0,5):
        jacobi(A, B, x0, x,eps)
        gass(A, B, x0, x, eps)
        print(y1data)
    A = np.array(y1data)
    B = np.array(y2data)
    index = np.arange(5)
    bar_width = 0.35
    plt.bar(index, A, bar_width,alpha=0.4, color='b')
    plt.bar(index + bar_width, B, bar_width,  alpha=0.4, color='r')
    x_labels = ['1', '2', '3', '4', '5']
    plt.xticks(index + bar_width / 2,x_labels)
    plt.show()
    #第二问
    y1data = []
    B1=[]
    for i in range(1, 21):
        B1.append(10)
    #print(B1)
    i=1
    while i<=8:
        data = np.array(C)
        row, col = np.diag_indices_from(data)
        data[row, col] = 3*i
        print(data)
        jacobi(data, B1, x0, x, eps)
        i+=1
    x_axis_data = [1, 2, 3, 4, 5,6,7,8]
    y_axis_data =np.array(y1data)
    plt.plot(x_axis_data, y_axis_data, 'ro-', color='#4169E1', alpha=0.8, linewidth=1)
    plt.xlabel('x轴数字')
    plt.ylabel('y轴数字')
    plt.show()

p115 上机四 4.1 拉格朗日插值：

def lagrange(n,x,y,xt,s=0):
    for i in range(n):
        p=1
        for j in range(i):
            p=p*(xt-x[j])/(x[i]-x[j])
        for j in range(i+1,n):
            p=p*(xt-x[j])/(x[i]-x[j])
        s=s+p*y[i]
    return s
if __name__=='__main__':
    x=[11,12,13]
    y=[0.190809,0.207912,0.224951]
    n=len(x)
    xt=11.5
    print(lagrange(n,x,y,xt))

p115 4.1 牛顿迭代法程序

def newton(x,y,n,t):
    for i in range(1,n+1):
        for j in range(n,i-1,-1):
            y[j]=(y[j]-y[j-1])/(x[j]-x[j-i])
    p,s=1,y[0]
    for i in range(1,n+1):
        p=p*(t-x[i-1])
        s=s+y[i]*p
    return s

if __name__=='__main__':
    x=[11,12,13]
    n=len(x)-1
    y=[0.190809,0.207912,0.224951]
    t=11.5
    answer=newton(x,y,n,t)
    print(answer)

p115上机四 4.2最小二拟合

import numpy as np
from scipy import linalg


def LU_decomposition(A, b):
    n = len(A[0])
    L = np.zeros([n, n])
    U = np.zeros([n, n])
    for i in range(n):
        L[i][i] = 1
        if i == 0:
            U[0][0] = A[0][0]
            for j in range(1, n):
                U[0][j] = A[0][j]
                L[j][0] = A[j][0] / U[0][0]
        else:
            for j in range(i, n):  # U
                temp = 0
                for k in range(0, i):
                    temp = temp + L[i][k] * U[k][j]
                U[i][j] = A[i][j] - temp
            for j in range(i + 1, n):  # L
                temp = 0
                for k in range(0, i):
                    temp = temp + L[j][k] * U[k][i]
                L[j][i] = (A[j][i] - temp) / U[i][i]

    Z = linalg.solve(L, b)
    x = linalg.solve(U, Z)
    return x


def nihe(x, y):
    # ans=np.zeros(2)
    xt = x.transpose()  # x1为x的转置
    x1 = np.matmul(xt, x)  # 矩阵相乘
    y1 = np.matmul(xt, y)
    # print(x1,y1)
    x = LU_decomposition(x1, y1)
    print(f'线性拟合曲线方程：y={x[0]}x+{x[1]}')


if __name__ == '__main__':
    x = np.array([[1, 1], [2, 1], [3, 1], [4, 1], [5, 1]])
    y = np.array([4, 4.5, 6, 8, 8.5])
    ans = np.zeros(2)
    nihe(x, y)

p145上机五 5.1复化梯形复化辛浦生

import math
import numpy as np

def f2(x):
    if (float(np.fabs(x)) < 1e-15):
        return 0
    y=math.sqrt(1-math.exp(-x))/x
    return y
def f3(x):
    y=math.sqrt(1-math.sin(x)**2)
    return y
def FHTx(f,a,b,n):
    ti=0.0
    h=(b-a)/n
    ti=f(a)+f(b)
    for k in range(1,int(n)):
        xk=a+k*h
        ti = ti + 2 * f(xk)
    FHTx = ti*h/2
    return FHTx
def FHXPs(f,a,b,n):
    si=0.0
    h = (b - a) / (2 * n)
    si=f(a)+f(b)
    for k in range(1,int(n)):
        xk = a + k * 2 * h
        si = si + 2 * f(xk)
    for k in range(int(n)):
        xk = a + (k * 2 + 1) * h
        si = si + 4 * f(xk)
    FHXPs = si*h/3
    return FHXPs
def main():
    T1=[]
    T2=[]
    a = math.pi/6
    b = 1
    a = float(a)
    b = float(b)
    i=0
    n=2
    n= float(n)
    while 1:
        T1.append(FHTx(f2,a,b,n))
        T1.append(FHTx(f2,a,b,n*2))
        n=n*2
        if(T1[i+1]-T1[i]<3*0.5*10**(-6)):
            print('复化梯形公式计算结果为',T1[i+1],'  n=',n)
            break
        i = i + 1
    j=0
    n=2
    n= float(n)
    while 1:
        T2.append(FHXPs(f2, a, b, n))
        T2.append(FHXPs(f2, a, b, n * 2))
        n = n * 2
        if (T2[j + 1] - T2[j] < 0.5 * 10 ** (-6)):
            print('复化辛浦生公式计算结果为',  T2[j + 1],'  n=',n)
            break
        j = j + 1


if __name__ == '__main__':
    main()

p145 上机五 5.3 高斯-勒让德求积

#高斯-勒让德求积公式
from sympy import *
from scipy.special import perm,comb  #排列,组合
import math
x,t = symbols("x,t")
#积分区间
a = 1
b = 3
#需要求积的目标函数
def f(x):
    f =math.exp(x)*math.sin(x)
    return f

n = 3  #n次多项式正交，n越大精度越高(n=0,1,2,...)
#勒让德多项式
def L(n):
    df = diff((x ** 2 - 1) ** (n + 1), x, n + 1)
    # Python内置阶乘函数factorial
    L = 1 /2**(n+1)/factorial(n+1) * df
    return L

#高斯点x求取
def Gauss_point(n):
    x_k_list = solve(L(n))   #求得零点解集
    return x_k_list

#求积系数A
def Quadrature_coefficient(x_k_list):
    A_list = []
    for x_k in x_k_list:
        A = 2/(1-x_k**2)/(diff(L(n),x,1).subs(x,x_k))**2
        A_list.append(A)
    return A_list

result = 0
x_k_list = Gauss_point(n)
A_list = Quadrature_coefficient(Gauss_point(n))
sum = len(A_list)
#区间变换
if a == -1 and b == 1:
    for i in range(sum):
        result += (A_list[i] * f(x_k_list[i])).evalf()
    print(result)
#将求求粉公式中的区间(a,b)转换为[-1,1]
else:
    for i in range(sum):
        X = (b-a)/2 * x_k_list[i] + (a+b)/2  #区间变换
        result += (b-a)/2 * (A_list[i] * f(X)).evalf()
    print(result)

p169 上机六 6.1 改进的欧拉法

def f1(x, y):
    res = x**2 + x**3 * y
    return res
def f2(x,y):
    res=y-2*x/y
    return res

def oula(f,xlist, ylist, h,a,b):
    n=(b-a)/h
    for i in range(int(n)):
        k1=h*f(xlist[i],ylist[i])
        k2=h*f(xlist[i]+h,ylist[i]+k1)
        y=ylist[i]+0.5*k1+0.5*k2
        ylist.append(y)
        xlist.append(xlist[i]+h)
if __name__=='__main__':
    xlists=[1]
    ylists=[1]
    h,a,b=0.1,1,1.5
    oula(f1,xlist=xlists,ylist=ylists,h=h,a=a,b=b)
    for i in range(len(xlists)):
        print(f'i={i},xi={xlists[i]},yi={ylists[i]}')