Untitled
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