one-file-projects/rk3.py

#!/usr/bin/env python
import math
import matplotlib.pyplot as plt

def rk3(h,tn,yn,ys,param=[]):
    k1 = ys(tn,yn,*param)
    k2 = ys(tn+h/2,yn+h*k1/2,*param)
    k3 = ys(tn+h,yn - h*k1 + 2*h*k2, *param)
    return tn+h, yn + h*(k1/3 + k2/3 + k3/3)

def main():
    U = lambda tn,f : math.sin(2*math.pi*tn*f)
    R = 5e3
    C = 1e-6
    us = lambda tn, un, f: (U(tn,f) - un)/(R*C)
    fc = 1.0/(2*math.pi*R*C)
    frequencies = [ 2**(i) * fc for i in range(-5,8)]
    ucmax_list = []
    ucmin_list = []
    for f in frequencies:
        print("Simulation f = %.3f Hz" % f)
        ucmax = 0
        ucmin = 0
        tn = 0
        un = 0
        for i in range(int(1000*f)):
            tn, un = rk3(1/(1000*f), tn, un, us, [f])
            ucmax = max(un,ucmax)
            ucmin = min(un,ucmin)
        ucmax_list.append(ucmax)
        ucmin_list.append(ucmin)

    plt.plot(frequencies, ucmax_list, 'r', frequencies, ucmin_list, 'b').show()


if __name__ == '__main__':
    main()