1 #!/bin/python3 2 # -*-coding:utf-8 -*- 3 import math 4 import numpy as np 5 6 #定义一个欧拉算法的类，从而实现不同步长的引用 7 class Euler: 8 y_list=[] #定义一个空列表来实现y值的存储 9 def __init__(self, h=0.1, y0=1,): #初始化Euler类的方法 10 self.h = h 11 self.y0 = y0 12 self.y = y0 13 self.n = 1/h 14 self.y_list = Euler.y_list 15 16 def euler(self, y0, h,): #定义Euler算法的过程 17 x = 0 18 y = y0 19 n = 1/h 20 for i in range(int(self.n+1)): 21 y_dere = 1 + math.log(x + 1) 22 y += self.h * y_dere 23 self.y_list = Euler.y_list.append("%.10f" % y) 24 25 # print("%.3f" % x, "%.10f" % y) 26 x += self.h 27 # print(Euler.y_list,) 28 return np.linspace(0, 1, n+1,dtype=float), Euler.y_list #为了引用结果列表，从而返回两个列表 29 30 31 #定义一个实现龙格库塔方法 32 class RungeKutta: 33 y_list = [] #同理，设置一个空列表来储存y值 34 def __init__(self, h=0.1, y0=1,): 35 self.h = h 36 self.y0 = y0 37 38 def make_ks(self, x_k,y_k): 39 40 def y_ders(x_k): 41 y_der = 1 + math.log(x_k + 1) 42 return y_der 43 #计算K1,K2,K3,K4 44 k1 = y_ders(x_k) 45 k2 = y_ders(x_k + self.h/2) 46 k3 = y_ders(x_k + (self.h/2)) 47 k4 = y_ders(x_k + self.h) 48 y_k = y_k+self.h/6*(k1+(2-math.sqrt(2))*k2 + (2 + math.sqrt(2))*k3 + k4) 49 return x_k, y_k 50 51 def make(self, h, y0): 52 n = int(1/self.h) 53 y_k = y0 54 x_k = 0 55 for i in range(n): #用循环遍历来实现y值的计算 56 x_k, y_k = self.make_ks(x_k,y_k) 57 x_k += h 58 # print("%.3f" % x_k, "%.10f" % y_k) 59 RungeKutta.y_list.append("%.10f" % y_k) 60 61 return np.linspace(0, 1, 1/h+1, dtype=float), RungeKutta.y_list 62 63 64 #下面将收集计算的值，首先输出Euler法，再输出rungkutta，x,y数据的顺序是：x,y 65 print("Euler", "h=0.1") 66 a1 = Euler() 67 ax1, ay1 = a1.euler(h=0.1, y0=1) 68 print(ax1, " ", ay1) 69 print("----------------------------------------------------------------------------") 70 print("Euler", "h=0.1") 71 print("h=0.01") 72 a2 = Euler() 73 ax2, ay2 = a2.euler(h=0.01, y0=1) 74 print(ax2, " ", ay2) 75 print("----------------------------------------------------------------------------") 76 print("Euler", "h=0.1") 77 print("h=0.2") 78 a3 = Euler() 79 ax3, ay3 = a3.euler(h=0.2, y0=1) 80 print(ax3, " ", ay3) 81 print("----------------------------------------------------------------------------") 82 83 print("Rungkutta", "h=0.2") 84 c1 = RungeKutta(h=0.2, y0=1) 85 cx1, cy1 = c1.make(h=0.2, y0=1) 86 print(cx1, " ", cy1) 87 print("-----------------------------------------------------------------------------") 88 89 print("Rungkutta", "h=0.1") 90 c2 = RungeKutta(h=0.1, y0=1) 91 cx2, cy2 = c2.make(h=0.1, y0=1) 92 print(cx2, " ", cy2) 93 print("-----------------------------------------------------------------------------") 94 95 print("Rungkutta", "h=0.01") 96 c3 = RungeKutta(h=0.01, y0=1) 97 cx3, cy3 = c3.make(h=0.01, y0=1) 98 print(cx3, " ", cy3) 99 print("-----------------------------------------------------------------------------")
以上的代码是实现算法的过程，以及输出的部分，还可以加入自定义的文件读取的办法，以及其他存储方式保存数据结果。
》》》》》》》》》》》》》》》》》》》》》》
程序2.0版本：
》》》》》》》》》》》》》》》》》》》》》》
# -*-coding:utf-8 -*- import math import numpy as np import matplotlib.pyplot as plt # 定义一个欧拉算法的类，从而实现不同步长的引用 class Euler: # y_list = [] # 定义一个空列表来实现y值的存储 def __init__(self, h, y0): # 初始化Euler类的方法 self.h = h self.y0 = y0 self.y = y0 self.n = 1 / self.h self.x = 0 def euler(self): # 定义Euler算法的过程 # x = 0 # y = self.y0 # n = 1 / self.h for i in range(int(self.n + 1)): y_dere = 1 + math.log(self.x + 1) self.y += self.h * y_dere self.y_list = Euler_y_list.append("%.10f" % self.y) self.x += self.h # print(Euler.y_list,) return np.linspace(0, 1, self.n + 1, dtype=float), Euler_y_list # 为了引用结果列表，从而返回两个列表 # 定义一个实现龙格库塔方法 class RungeKutta: def __init__(self, h, y0=1): self.h = h self.y0 = y0 def make_ks(self, x_k, y_k): def y_ders(x_k): y_der = 1 + math.log(x_k + 1) return y_der # 计算K1,K2,K3,K4 k1 = y_ders(x_k) k2 = y_ders(x_k + self.h / 2) k3 = y_ders(x_k + (self.h / 2)) k4 = y_ders(x_k + self.h) y_k = y_k + self.h / 6 * (k1 + (2 - math.sqrt(2)) * k2 + (2 + math.sqrt(2)) * k3 + k4) return x_k, y_k def make(self): n = int(1 / self.h) y_k = self.y0 x_k = 0 for i in range(n): # 用循环遍历来实现y值的计算 x_k, y_k = self.make_ks(x_k, y_k) x_k += self.h RungeKutta_y_list.append(y_k) return np.linspace(0, 1, n+1, endpoint=True,dtype=float), RungeKutta_y_list class Plots: def __init__(self, x, y, title): self.x = list(x) self.y = list(y) self.title = title def plots(self): plt.figure(figsize=(8, 4)) plt.scatter(self.x, self.y, label="y'=1+ln(x+1), y0=1 (0<x<=1) ", color="red", linewidth=2) plt.xlabel("x") plt.ylabel("y") plt.title(self.title) # plt.ylim(min(self.y), max(self.y)) plt.ylim() plt.legend() plt.show() if __name__ == "__main__": while True: name = input("本程序提供Euler算法和Rungekutta算法 选择算法时输入exit退出 请输入需要的算法名：") try: a = float(input("请输入所需要的h值：")) if a =="exit" or a=="Exit": break except ValueError as ve: print(ve) if name == "Euler": Euler_y_list =[] l1 = Euler(h=a, y0=1) l1x, l1y = l1.euler() with open("Euler.txt","a+") as f: f.write("h=%s " % l1.h) f.write("x值--------------y值 ") for i in range(len(l1x)): print("x是：%.5f y是：%.10f" % (float(l1x[i]),float(l1y[i]))) a,b = str(l1x[i]),str(l1y[i]) f.write(a+" "*6+b+" ") l1p = Plots(l1x, l1y, name) l1p.plots() print("请储存好你的x值，和y值") del l1, l1p, l1x, l1y Euler.y_list = [] elif name == "Rungekutta": RungeKutta_y_list = [0] l2 = RungeKutta(h=a) l2x, l2y = l2.make() with open("Rungekutta.txt", "a+") as f: f.write("h=%s " % l2.h) f.write("x值--------------y值 ") for i in range(len(l2x)): print("x是：%.5f y是：%.10f" % (l2x[i], l2y[i])) a, b = str(l2x[i]), str(l2y[i]) f.write(a + " " * 6 + b + " ") l2p = Plots(l2x, l2y, name) l2p.plots() print("请储存好你的x值，和y值") del l2, l2p elif name =="exit": break else: print("请输入正确的算法名!!!!")
加入了退出和写入到文件的功能