import numpy as np import matplotlib.pyplot as plt ''' Implicit vs explicit Euler method for stiff problems function oppgave_6_6_2 ''' def ydot(t,y): # right-hand side of differential equation #z = z=6*y-3*y**2 # oppgave (a) z = 10*y**3-10*y**4 # oppgave (b) return z def ydot_prime(t,y): # derivative of the right-hand side with respect to y # (used for implicit methods) #z = 6-6*y# oppgave 6.6.2 (a) z = 30*y**2-40*y**3 # oppgave 6.6.2 (b) 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 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): 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 oppgave_6_6_2(y_0,a,b,n): h = (b-a)/n t = a ye = np.zeros(n+1) yi = np.zeros(n+1) # start the integration loop ye[0] = y_0 yi[0] = y_0 for i in range(n): ye[i+1] = eulerstep(t,ye[i],h) yi[i+1] = impeulerstep(t,yi[i],h) t = t+h # plot the explicit vs implicit approximation plt.figure(0) plt.clf() fs = 24 #font size plt.plot(np.arange(a,b+h,h),yi,'b-',lw=3,markersize=10) plt.plot(np.arange(a,b+h,h),ye,'r-',lw=3,markersize=10) plt.xlabel('t',fontsize=fs) plt.ylabel('y(t)',fontsize=fs) plt.legend(['Explicit Euler','Implicit Euler']) plt.title('Explicit vs implicit Euler, h=%e'%(h)) plt.show() if __name__=='__main__': y_0 = 0.5 a = 0. b = 20. n = 66 oppgave_6_6_2(y_0,a,b,n)