Lindstedt Perturbation Method

(t, epsilon) = var('t, epsilon') 

# Regular perturbation approximation # Divided by epsilon thetaR1(t) = cos(t) thetaR3(t) = thetaR1(t) + epsilon^2 / 192 * (cos(t) - cos(3 * t) + 12 * t * sin(t)) # Display the above html(r'1-Term Regular Approx: $\theta_{R1}(t) = %s$' % latex(epsilon * thetaR1(t))) html(r'2-Term Regular Approx: $\theta_{R3}(t) = %s$' % latex(epsilon * thetaR3(t))) 

1-Term Regular Approx: \theta_{R1}(t) = \epsilon \cos\left(t\right)
2-Term Regular Approx: \theta_{R3}(t) = \frac{1}{192} \, {\left({\left(12 \, t \sin\left(t\right) - \cos\left(3 \, t\right) + \cos\left(t\right)\right)} \epsilon^{2} + 192 \, \cos\left(t\right)\right)} \epsilon

1-Term Regular Approx: \theta_{R1}(t) = \epsilon \cos\left(t\right)
2-Term Regular Approx: \theta_{R3}(t) = \frac{1}{192} \, {\left({\left(12 \, t \sin\left(t\right) - \cos\left(3 \, t\right) + \cos\left(t\right)\right)} \epsilon^{2} + 192 \, \cos\left(t\right)\right)} \epsilon

# Lindstedt approximation tau0(t) = t tau2(t) = (1 - epsilon^2 / 16) * t # Divided by epsilon thetaL1(tau) = cos(tau) thetaL3(tau) = thetaL1(tau) + epsilon^2 / 192 * cos(tau) # Display the above html(r'1-Term Strained Time: $\tau_0 = %s$' % latex(tau0(t))) html(r'2-Term Strained Time: $\tau_2 = %s$' % latex(tau2(t))) html(r'1,1-Term Lindstedt Approx: $\theta_{L1}(\tau_0) = %s$' % latex(epsilon * thetaL1(tau0(t)))) html(r'1,2-Term Lindstedt Approx: $\theta_{L1}(\tau_2) = %s$' % latex(epsilon * thetaL1(tau2(t)))) html(r'2,2-Term Lindstedt Approx: $\theta_{L3}(\tau_2) = %s$' % latex(epsilon * thetaL3(tau2(t)))) 

1-Term Strained Time: \tau_0 = t
2-Term Strained Time: \tau_2 = -\frac{1}{16} \, {\left(\epsilon^{2} - 16\right)} t
1,1-Term Lindstedt Approx: \theta_{L1}(\tau_0) = \epsilon \cos\left(t\right)
1,2-Term Lindstedt Approx: \theta_{L1}(\tau_2) = \epsilon \cos\left(-\frac{1}{16} \, {\left(\epsilon^{2} - 16\right)} t\right)
2,2-Term Lindstedt Approx: \theta_{L3}(\tau_2) = \frac{1}{192} \, {\left(\epsilon^{2} \cos\left(-\frac{1}{16} \, {\left(\epsilon^{2} - 16\right)} t\right) + 192 \, \cos\left(-\frac{1}{16} \, {\left(\epsilon^{2} - 16\right)} t\right)\right)} \epsilon

1-Term Strained Time: \tau_0 = t
2-Term Strained Time: \tau_2 = -\frac{1}{16} \, {\left(\epsilon^{2} - 16\right)} t
1,1-Term Lindstedt Approx: \theta_{L1}(\tau_0) = \epsilon \cos\left(t\right)
1,2-Term Lindstedt Approx: \theta_{L1}(\tau_2) = \epsilon \cos\left(-\frac{1}{16} \, {\left(\epsilon^{2} - 16\right)} t\right)
2,2-Term Lindstedt Approx: \theta_{L3}(\tau_2) = \frac{1}{192} \, {\left(\epsilon^{2} \cos\left(-\frac{1}{16} \, {\left(\epsilon^{2} - 16\right)} t\right) + 192 \, \cos\left(-\frac{1}{16} \, {\left(\epsilon^{2} - 16\right)} t\right)\right)} \epsilon

# ODE for numerical solution def RHS(t, Y): (theta, thetaDot) = Y dTheta = thetaDot dThetaDot = - sin(theta) return (dTheta, dThetaDot) # Initial condition Y0(epsilon) = [epsilon, 0] # Set up solver T = ode_solver(function = RHS, algorithm = 'rk8pd') 

# Interactively plot solution & approximations # Set up controls for interactive plots def showCheckBox(*args, **kwargs): return checkbox(*args, **kwargs) def colorSelector(*args, **kwargs): return color_selector(*args, label = 'Color', hide_box = True, widget = 'jpicker', **kwargs) # The same transparency for each curve transparency = 0.5 lineStyles = ('-', '--', '-.', ':') def lineStyleSelector(*args, **kwargs): return selector(*args, values = lineStyles, label = 'Style', **kwargs) thicknesses = (1, 2, 3) def thicknessSelector(*args, **kwargs): return selector(*args, values = thicknesses, label = 'Thickness', **kwargs) # Make the plots @interact(layout = [['log10Epsilon'], ['tRange'], ['NShown', 'NColor', 'NStyle', 'NThickness'], ['R1Shown', 'R1Color', 'R1Style', 'R1Thickness'], ['R3Shown', 'R3Color', 'R3Style', 'R3Thickness'], ['L10Shown', 'L10Color', 'L10Style', 'L10Thickness'], ['L12Shown', 'L12Color', 'L12Style', 'L12Thickness'], ['L32Shown', 'L32Color', 'L32Style', 'L32Thickness']]) def plotSeries(log10Epsilon = slider(-3, 0, default = log(1. / 3, 10), label = '$\log_{10} \epsilon$'), tRange = range_slider(vmin = 0, vmax = 1000, default = (0, 200), label = '$t$'), NShown = showCheckBox(label = 'Numerical Solution'), NColor = colorSelector(default = (0, 0, 0)), NStyle = lineStyleSelector(default = '--'), NThickness = thicknessSelector(), R1Shown = showCheckBox(label = '1-Term Regular Approx'), R1Color = colorSelector(default = (1, 0, 0)), R1Style = lineStyleSelector(), R1Thickness = thicknessSelector(), R3Shown = showCheckBox(label = '2-Term Regular Approx'), R3Color = colorSelector(default = (1, 1, 0)), R3Style = lineStyleSelector(), R3Thickness = thicknessSelector(), L10Shown = showCheckBox(default = False, label = '1,1-Term Lindstedt Approx'), L10Color = colorSelector(default = (0, 1, 0)), L10Style = lineStyleSelector(), L10Thickness = thicknessSelector(), L12Shown = showCheckBox(label = '1,2-Term Lindstedt Approx'), L12Color = colorSelector(default = (0, 1, 1)), L12Style = lineStyleSelector(), L12Thickness = thicknessSelector(), L32Shown = showCheckBox(default = False, label = '2,2-Term Lindstedt Approx'), L32Color = colorSelector(default = (1, 0, 1)), L32Style = lineStyleSelector(), L32Thickness = thicknessSelector()): epsilon = 10^log10Epsilon html(r'$$\epsilon = %s$$' % latex(epsilon)) (tMin, tMax) = tRange nPoints = int(100 / 2 / pi * (tMax - tMin)) t_span = srange(tMin, tMax, (tMax - tMin) / (nPoints - 1), include_endpoint = True) if tMin > 0: t_span = [0] + t_span T.ode_solve(y_0 = Y0(epsilon), t_span = t_span) tN = [s[0] for s in T.solution if s[0] >= tMin] thetaN = [s[1][0] / epsilon for s in T.solution if s[0] >= tMin] solutionPlots = [] errorPlots = [] if R1Shown: solutionPlots.append(plot(thetaR1(t)(epsilon = epsilon), t, tMin, tMax, plot_points = nPoints, color = R1Color, alpha = transparency, linestyle = R1Style, thickness = R1Thickness)) err = [(thetaN[i] - thetaR1(tN[i])(epsilon = epsilon)) for i in range(len(tN))] errorPlots.append(line(zip(tN, err), color = R1Color, alpha = transparency, linestyle = R1Style, thickness = R1Thickness)) if R3Shown: solutionPlots.append(plot(thetaR3(t)(epsilon = epsilon), t, tMin, tMax, plot_points = nPoints, color = R3Color, alpha = transparency, linestyle = R3Style, thickness = R3Thickness)) err = [(thetaN[i] - thetaR3(tN[i])(epsilon = epsilon)) for i in range(len(tN))] errorPlots.append(line(zip(tN, err), color = R3Color, alpha = transparency, linestyle = R3Style, thickness = R3Thickness)) if L10Shown: solutionPlots.append(plot(thetaL1(tau0(t))(epsilon = epsilon), t, tMin, tMax, plot_points = nPoints, color = L10Color, alpha = transparency, linestyle = L10Style, thickness = L10Thickness)) err = [(thetaN[i] - thetaL10(tau0(tN[i]))(epsilon = epsilon)) for i in range(len(tN))] errorPlots.append(line(zip(tN, err), color = L10Color, alpha = transparency, linestyle = L10Style, thickness = L10Thickness)) if L12Shown: solutionPlots.append(plot(thetaL1(tau2(t))(epsilon = epsilon), t, tMin, tMax, plot_points = nPoints, color = L12Color, alpha = transparency, linestyle = L12Style, thickness = L12Thickness)) err = [(thetaN[i] - thetaL1(tau2(tN[i]))(epsilon = epsilon)) for i in range(len(tN))] errorPlots.append(line(zip(tN, err), color = L12Color, alpha = transparency, linestyle = L12Style, thickness = L12Thickness)) if L32Shown: solutionPlots.append(plot(thetaL3(tau2(t))(epsilon = epsilon), t, tMin, tMax, plot_points = nPoints, color = L32Color, alpha = transparency, linestyle = L32Style, thickness = L32Thickness)) err = [(thetaN[i] - thetaL3(tau2(tN[i]))(epsilon = epsilon)) for i in range(len(tN))] errorPlots.append(line(zip(tN, err), color = L32Color, alpha = transparency, linestyle = L32Style, thickness = L32Thickness)) if NShown: solutionPlots.append(line(zip(tN, thetaN), color = NColor, alpha = transparency, linestyle = NStyle, thickness = NThickness)) show(sum(solutionPlots), axes_labels = (r'$t$', r'$\theta / \epsilon$')) show(sum(errorPlots), axes_labels = (r'$t$', 'Abs. Err.'))