Untitled

unknown
python
2 months ago
3.3 kB
2
Indexable
Never
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()