Solving homogeneous mass-spring systems

%hide %html An interact is a small applet that allows you to change the parameters and see results. This interact solves the differential equation modeling a spring system which has no forcing function. 

%hide #auto RF=RealField(10) var("t") t0=float(0) 

CPU time: 0.00 s,  Wall time: 0.00 s

CPU time: 0.00 s,  Wall time: 0.00 s

%hide %cython #auto #from math import exp cdef extern from "math.h": double exp(double x) double cos(double x) double sin(double x) double atan2(double y, double x) cpdef list solve_diffeq(double a, double b, double c, double t0, double y0, double y0p): """ Returns a list starting with discriminant if discriminant > 0, then the list is [disc, r1, r2, c1, c2] if discriminant == 0, then the list is [disc,r1,c1,c2] if discriminant < 0, then the list is [disc, lambda, mu, c1, c2,r,delta] """ cdef list result = [] cdef double discriminant = b*b-4*a*c cdef double const[2] result.append(discriminant) cdef double r1, r2, pos_disc_root, y1t0, y1difft0, y2t0, y2difft0, denominator, c1, c2 cdef double lam, mu, elam, sinmu, cosmu, r, delta if discriminant > 0: pos_disc_root = sqrt(discriminant) r1 = (-b + pos_disc_root)/(2*a) r2 = (-b - pos_disc_root)/(2*a) result.append(r1) result.append(r2) y1t0 = exp(r1*t0) y1difft0 = r1*y1t0 y2t0 = exp(r2*t0) y2difft0 = r2*y2t0 denominator = y1t0*y2difft0 - y2t0*y1difft0 c1 = (y0*y2difft0 - y0p*y2t0)/denominator c2 = (-y0*y1difft0 + y0p*y1t0)/denominator result.append(c1) result.append(c2) elif discriminant == 0: r1 = -b/(2*a) result.append(r1) y1t0 = exp(r1*t0) y1difft0 = r1*exp(r1*t0) y2t0 = t0*y1t0 y2difft0 = y1t0 + r1*y2t0 denominator = y1t0*y2difft0 - y2t0*y1difft0 c1 = (y0*y2difft0 - y0p*y2t0)/denominator c2 = (-y0*y1difft0 + y0p*y1t0)/denominator result.append(c1) result.append(c2) else: lam = -b/(2*a) mu = sqrt(-discriminant)/(2*a) result.append(lam) result.append(mu) elam = exp(lam*t0) cosmu = cos(mu*t0) sinmu = sin(mu*t0) y1t0 = elam*cosmu y1difft0 = elam*(lam*cosmu-mu*sinmu) y2t0 = elam*sinmu y2difft0 = elam*(lam*sinmu+mu*cosmu) denominator = y1t0*y2difft0 - y2t0*y1difft0 c1 = (y0*y2difft0 - y0p*y2t0)/denominator c2 = (-y0*y1difft0 + y0p*y1t0)/denominator result.append(c1) result.append(c2) r = sqrt(c1*c1+c2*c2) delta = atan2(c2,c1) result.append(r) result.append(delta) return result 

CPU time: 0.04 s,  Wall time: 4.25 s

CPU time: 0.04 s,  Wall time: 4.25 s

%hide #auto @interact def _(initial=input_grid(1,3,default=[-1/12,2,5],label="y(0), y'(0),t_max",width=10), abc=input_grid(1,3,default=[3/32,0,12],label="m,gamma,k", width=10)): y0,y0p,tmax = initial[0] a,b,c = abc[0] plots = [] a = float(a) b = float(b) c = float(c) soln_nums = solve_diffeq(a,b,c,t0,y0,y0p) discriminant = soln_nums[0] if discriminant > 0: r1,r2,c1,c2 = soln_nums[1:] soln = c1*e^(r1*t)+c2*e^(r2*t) linear_combo = ["%se^{%st}"%(RF(constant),RF(rval)) for constant,rval in ((c1,r1),(c2,r2)) if constant !=0] soln_n = "$y(t) = " + '+'.join(linear_combo) + "$" elif discriminant == 0: r1,c1,c2 = soln_nums[1:] soln = (c1+c2*t)*e^(r1*t) linear_combo = ["%s%s"%(RF(constant),y) for constant,y in ((c1,"e^{%st}"%RF(r1)),(c2,"te^{%st}"%RF(r1))) if constant !=0] soln_n = "$y(t) = " + '+'.join(linear_combo) + "$" else: lam,mu,c1,c2,R,delta = soln_nums[1:] soln = e^(lam*t)*(c1*cos(mu*t)+c2*sin(mu*t)) linear_combo = ["%s%s"%(RF(constant),y) for constant,y in ((c1,"\cos(%st)"%RF(mu)),(c2,"\sin(%st)"%RF(mu))) if constant !=0] soln_n = "$y(t) = e^{%st}("%RF(lam) + '+'.join(linear_combo) + ")$" soln2 = R*e^(lam*t)*cos(mu*t-delta) soln2_n = "%se^{%st}\cos(%st-%s)"%(RF(R),RF(lam),RF(mu),RF(delta)) plots.extend([plot(R*e^(lam*t),(t,0,tmax),rgbcolor=(0.1,0.1,0.1)), plot(-R*e^(lam*t),(t,0,tmax),rgbcolor=(0.1,0.1,0.1))]) line1 = "$\\left(%s \\right)y''+\\left(%s \\right)y'+\\left(%s \\right)y=0; \quad$"%(latex(a),latex(b),latex(c)) line1 += "$y(%s)=%s,\quad y'(%s)=%s$"%(latex(t0), latex(y0), latex(t0), latex(y0p)) line1 += "$\quad$ discriminant = $%s$"%latex(discriminant) line2 = soln_n if discriminant < 0: line2 += "$\quad = \quad "+soln2_n+"$" html(line1) html(line2) plots.append(plot(soln,(t,0,tmax))) show(sum(plots))