Untitled

 avatar
unknown
python
2 years ago
3.3 kB
6
Indexable
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
import math
import time

print("введите количество разбиений N = ", end="")
N = int(input())
print("введите параметр p = ", end = "")
p = float(input())
x = []
a = 1.
b = 2.
h = 1./N
eps = 0.000001
for i in range(N):
    x.append(1+i*h)


#def fi(x):
#    return np.sin(np.pi*x)
#def d_fi(x):
#    return np.pi*np.cos(np.pi*x)

#def fi(x):
#    return np.log((2*x)+1)
#def d_fi(x):
#    return 2/((2*x)+1)

def fi(x):
    return np.log(x)
def d_fi(x):
    return 1/(x)

def F(x, u):
    return d_fi(x) + p*(u-fi(x))**3
def dF(x, u):
    return p*3*((u-fi(x))**2)

def H(x, u_old, u_new):
    return u_new - u_old - h*F(x,u_new)

def dH(x, u):
    return 1-h*dF(x,u)

def Nuton(x,u):
    u_old = (a+b)/2
    u_new = u_old
    k=0
    while((abs(H(x, u, u_new)))>eps):
        if(dH(x, u_new)==0):
            u_old = u_old +eps
        else:
            k += 1
            u_old = u_new
            u_new = u_old -H(x, u, u_old)/dH(x, u_old)

    return u_new

u = []
u.append(0)


#Неявный метод Эйлера
start = time.time()
for i in range(1,N):
    u.append(Nuton(x[i-1], u[i-1]))

Fun = []
for i in range(N):
    Fun.append(fi(x[i]))
plt.scatter(x, Fun,color='black')

end = time.time()
Eror = 0
#Вычисляем порядок точности для метода Эйлера
for i in range(N):
    e = abs(Fun[i] - u[i])
    if(Eror < e):
        Eror = e

plt.scatter(x, u,color='blue')
print("Количество узлов N =",N, "Ошибка Эйлер =",Eror,"time Эйлер",(end-start), "ms")

#Рунге-кутты четвертого порядка
start = time.time()
for j in range(1000):
    for i in range(N-1):
        k1 = h*F(x[i],u[i])
        k2 = h*F(x[i]+h/2,u[i]+k1/2)
        k3 = h*F(x[i]+h/2,u[i]+k2/2)
        k4 = h*F(x[i]+h,u[i]+k3)
        u[i+1] = u[i] + (k1+2*k2+2*k3+k4)/6
end = time.time()

Eror = 0
#Вычисляем порядок точности для метода Рунге-Кутты
for i in range(N):
    e = abs(Fun[i] - u[i])
    if(Eror < e):
        Eror = e



plt.scatter(x, u, color='green')
print("Количество узлов N =",N, "Ошибка Runge =",Eror,"time Runge",(end-start), "ms")

#Метод Аддамса четвертого порядка

start_time = time.time()
for j in range(1000):
    a1, a2, a3, a4 = F(x[0], u[0]), F(x[1], u[1]), F(x[2], u[2]), F(x[3], u[3])
    for i in range(3, N-1):
        u[i+1] = u[i] + (h/24)*(55*a4 - 59*a3 + 37*a2 - 9*a1)
        a1, a2, a3 = a2, a3, a4
        a4 = F(x[i+1], u[i+1])
end_time = time.time()
plt.scatter(x,u,color='red')

Eror = 0
#Вычисляем порядок точности для Аддамса
for i in range(N):
    e = abs(Fun[i] - u[i])
    if(Eror < e):
        Eror = e
print("Количество узлов N =",N, "Ошибка Аддамса =",Eror,"time Addams",(end_time-start_time), "ms")


plt.text(1, 2.5, "Эйлер - Blue")
plt.text(1, 2.0, "Рунге-кутты - Green")
plt.text(1, 1.5, "Adams - Red")
plt.text(1, 1.0, "решение - Black")

plt.show()
Editor is loading...
Leave a Comment