import numpy as np import matplotlib.pyplot as plt import scipy.special ''' Solve the stiff ODE y' = y^2 - y^3 using various explicit/implicit methods ''' def ydot(t,y): # right-hand side of differential equation z = y**2-y**3 return z def ydot_prime(t,y): # derivative of the right-hand side with respect to y # (used for implicit methods) z=2*y-3*y**2 return z def eulerstep(t,y,h): # one step of Euler's Method # Input: current time t, current value y, stepsize h # Output: approximate solution value at time t+h y=y+h*ydot(t,y) return y def trapstep(t,y,h): # one step of explicit Trapezoid Method # Input: current time t, current value y, stepsize h # Output: approximate solution value at time t+h f1 = ydot(t,y) f2 = ydot(t+h,y+h*f1) y = y + h/2*(f1+f2) return y def midstep(t,y,h): #one step of Midpoint Method #Input: current time t, current value y, stepsize h # Output: approximate solution value at time t+h k1= ydot(t,y) k2= ydot(t+h/2,y+h/2*k1) y = y + h*k2 return y def impeulerstep(t,y,h): # solve the equation y1-h*f(t+h,y1)=y # use Newton's method tol = 1.0e-10 MaxIt = 100 y1 = y for i in range(MaxIt): resid = y1 - h*ydot(t+h,y1) - y if ((abs(resid)/max(abs(y),1)) < tol): # success: stop iterating and return the solution return y1 # end if jac = 1 - h*ydot_prime(t+h,y1) y1 = y1 - resid/jac # end for print('Warning: Newton''s iteration did not converge at t=%e\n' % (t)) return y1 def imptrapstep(t,y,h): # solve the equation y1-h/2*f(t+h,y1)=y+h/2f(t,y) # use Newton's method tol = 1.0e-10 MaxIt = 100 y1 = y for i in range(MaxIt): resid = y1 - h/2*ydot(t+h,y1) - y - h/2*ydot(t,y) if ((abs(resid)/max(abs(y),1)) < tol): # success: stop iterating and return the solution return y1 #end if jac = 1 - h/2*ydot_prime(t+h,y1) y1 = y1 - resid/jac #end for print('Warning: Newton''s iteration did not converge at t=%e\n' % (t)) return y1 def int_flame(T,n,y_0): ''' Solve IVP: y' = y^2-y^3 y(0) = y_0 ''' h = T/n t = 0 y = np.zeros(n+1)+0j # start the integration loop y[0] = y_0 for i in range(n): y[i+1] = impeulerstep(t,y[i],h) t = t+h # exact solution alpha = 1/y_0 - 1 y_exact = lambda t: 1/(np.real(scipy.special.lambertw(alpha*np.exp(alpha-t)))+1) # plot the explicit vs exact solution plt.figure(0) plt.clf() fs = 24 #font size plt.plot(np.arange(0,T+h,h),y,'b-',lw=3,markersize=10) plt.plot(np.arange(0,T+h,h),y_exact(np.arange(0,T+h,h)),'r-',lw=3,markersize=10) plt.xlabel('t',fontsize=fs) plt.ylabel('y(t)',fontsize=fs) plt.figure(1) plt.clf() plt.plot(np.arange(0,T+h,h),y-y_exact(np.arange(0,T+h,h)),'rx-',lw=3,markersize=10) plt.xlabel('t',fontsize=fs) plt.ylabel('Error',fontsize=fs) plt.show() if __name__=='__main__': # smaller y_0 => stiffer solution y_0 = 1.0e-02 T = 2/y_0 n = 100 int_flame(T,n,y_0)