|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 |
import os import numpy as np import matplotlib.pyplot as plt os.chdir(r'D:\projects\wordpress\ex46') os.getcwd() #2nd-order Runge-Kutta methods with A=1/3 (type C) # dy/dx=exp(-2x)-2y # y(0)=0.1, interval x=[0,2], step size = h=0.2 def feval(funcName, *args): return eval(funcName)(*args) def RK2C(func, yinit, x_range, h): m = len(yinit) n = int((x_range[-1] - x_range[0])/h) x = x_range[0] y = yinit # Containers for solutions xsol = np.empty(0) xsol = np.append(xsol, x) ysol = np.empty(0) ysol = np.append(ysol, y) for i in range(n): k1 = feval(func, x, y) ypredictor = y + k1 * (3*h/4) k2 = feval(func, x+3*h/4, ypredictor) for j in range(m): y[j] = y[j] + (h/3)*(k1[j] + 2*k2[j]) x = x + h xsol = np.append(xsol, x) for r in range(len(y)): ysol = np.append(ysol, y[r]) return [xsol, ysol] def myFunc(x, y): dy = np.zeros((len(y))) dy[0] = np.exp(-2 * x) - 2 * y[0] return dy # ----------------------- h = 0.2 x = np.array([0, 2]) yinit = np.array([1.0/10]) [ts, ys] = RK2C('myFunc', yinit, x, h) dt = int((x[-1]-x[0])/h) t = [x[0]+i*h for i in range(dt+1)] yexact = [] for i in range(dt+1): ye = (1.0/10)*np.exp(-2*t[i]) + t[i]*np.exp(-2*t[i]) yexact.append(ye) plt.plot(ts, ys, 'r') plt.plot(t, yexact, 'b') plt.xlim(x[0], x[1]) plt.legend(["2nd Order RK, A=1/3 ", "Exact solution"], loc=1) plt.xlabel('x', fontsize=17) plt.ylabel('y', fontsize=17) plt.tight_layout() plt.savefig("example46.png", dpi=100) plt.show() plt.close() |
Recent Comments