# -*- coding: utf-8 -*- # 3.0 # # Numerical Solutions of ODE # # Euler and Implicit Euler # # The population model again # # In the previous ipython notebook, we solved an ordinary differential equation for population model: # %matplotlib inline from sympy import * import sympy as sp from pylab import * # define an ODE t = symbols('t') p = symbols('p',cls=Function) p = Function('p')(t) logisticODE = Eq(p.diff(t),0.03*p-0.0005*p**2) pprint(logisticODE) #solve an ODE #logisticSol = dsolve(logisticODE,p) #pprint(logisticSol) # substitude C1, so that p(0)=2 #C1 = Symbol('C1') #logisticFun = logisticSol.rhs.subs([(C1,Rational(29,60))]) # newer sympy versions fail to solve the system, older ones were working fine! # analytical solution logisticFun = 60/(1+29*functions.exp(-3*t/100)) pprint(logisticFun) # def plotDirectionField(X, Y, U, V, pname): """ simple plotting of direction field.""" M = np.sqrt(pow(U, 2) + pow(V, 2)) U = np.divide(U,M) V = np.divide(V,M) Q = quiver( X, Y, U, V, pivot='mid') title(pname) # plot the solution lam_logistic = lambdify((t), logisticFun) #set the limits x_min = 0. x_max = 300. y_min = 0. y_max = 60. Nt = 100 T = np.linspace(0.,x_max,Nt) P = [lam_logistic(Te) for Te in T] figure(figsize=(10, 10),dpi=120) plot(T, P, lw=2) # generate the direction field X,Y = meshgrid( arange(x_min-15.,x_max+15.,15.),arange(y_min-3.,y_max+3.,3.) ) U = np.copy(X) # adjust the direction field direction to the plot axis U.fill((y_max-y_min)/(x_max-x_min)) V = (0.03-0.0005*Y)*Y plotDirectionField(X, Y, U, V, "Logistic Growth") xlim([x_min,x_max]) ylim([y_min,y_max]) xlabel('t') ylabel('p') show() # # Preparation for the numerical solution # # Although we get an exact solution for this problem, we will solve it again numerically (so we # can compare the resulting approximation with the exact solution). # # First, for some integer number n, we divide the Interval [0,Tend] into n subintervals - here, we choose all subintervals of the same length delta_t = Tend/n: # TEnd = 100 n = 5 tau = TEnd/5 #initial condition p0 = 2 # # The numerical integration will compute values $pnum_k$, $k=1 \dotsc n+1$ for $p(tp_k)$ with $tp_k = (k-1) \tau$. We store the grid points $tp$ into a vector: # tp = [(k-1)*tau for k in range(1,n+2)] pprint(tp) # # We can draw the intervals with respect to t in the direction field: # # set the limits x_min = 0. x_max = 100. y_min = 0. y_max = 25. f1 = figure(figsize=(12, 12),dpi=120) # plot the analytical solution plot(T, P, lw=2) # plot the direction field X,Y = meshgrid( arange(x_min,x_max,5.),arange(y_min,y_max,1.) ) U = np.copy(X) # adjust the direction field direction to the plot axis U.fill((y_max-y_min)/(x_max-x_min)) V = (0.03-0.0005*Y)*Y plotDirectionField(X, Y, U, V, "Logistic Growth") # plot the intervals vlines = [[[i,i]]+[[y_min,y_max]] for i in tp] vlines = [j for i in vlines for j in i] plot(*vlines,c='r') xlim([x_min,x_max]) ylim([y_min,y_max]) xlabel('t') ylabel('p') show() # # Create an vector with $n+1$ elements $pp[1] \dotsc pp[n+1]$ (the method itself will not require this storage, but we save all intermediate results, so that we can plot them afterwards): # pp1 = zeros(n+1) pprint(pp1) # # We know the first entry of p from the initial condition: # pp1[0] = p0 # # Euler's method # # Explicit Euler # # To compute the other values $p_k, k=1 \dotsc n$, we use Euler's method, that approximates $p_{k+1}$ by $p_k + \tau f(t_k,p_k)$, where $f(t,p)$ is the right hand side of the differential equation: # # get the right hand side function from the logistic ODE # substitude the function p(t) by some symbol x and lambdify the result x = Symbol('x') lam_logistic_rhs = lambdify(x,logisticODE.rhs.subs(p,x)) for k in range(0,n): pp1[k+1] = pp1[k] + tau*lam_logistic_rhs(pp1[k]) pprint(pp1) f2 = figure(figsize=(12, 12),dpi=120) # plot analytical solution plot(T, P, lw=2,c='b') # plot numerical solution (Explicit Euler) plot(tp,pp1,lw=2,c='g') plotDirectionField(X, Y, U, V, "Logistic Growth") plot(*vlines,c='r') xlim([x_min,x_max]) ylim([y_min,y_max]) xlabel('t') ylabel('p') show() # # Implicit Euler # # Using the implicit Euler method, we no longer have an explicit formula for the solution in the next time step. # In principle, we would like to do something like this: # pp2 = zeros(n+1) pp2[0] = p0 for k in range(0,n): # note the change: pp[k+1] instead of pp[k] pp2[k+1] = pp2[k] + tau*lam_logistic_rhs(pp2[k+1]) pprint(pp2) # # However, this doesn't work - we have to solve an equation in every step, instead: # pp2 = zeros(n+1) pp2[0] = p0 logisticRHS = logisticODE.rhs.subs(p,x) pprint(logisticRHS) for k in range(0,n): # use pnew as the value in the next step # ft := subs(p(t)=pnew,rhs(popeq)); stepeq = Eq(x, Float(pp2[k]) + Float(tau)*logisticRHS) pprint(stepeq) stepeqSOL = sp.solve(stepeq,x) pprint(stepeqSOL) # the quadratic equation stepeq has two solutions # in this example, we can alsways take the bigger(positive) one pp2[k+1] = max(stepeqSOL) pprint(pp2) f2 = figure(figsize=(12, 12),dpi=120) # plot analytical solution plot(T, P, lw=2,c='b') # plot numerical solution (Explicit Euler) plot(tp,pp1,lw=2,c='g') # plot numerical solution (Implicit Euler) plot(tp,pp2,lw=2,c='y') plotDirectionField(X, Y, U, V, "Logistic Growth") plot(*vlines,c='r') xlim([x_min,x_max]) ylim([y_min,y_max]) xlabel('t') ylabel('p') show() # # Possible Problems # # Stability # # define an ODE stabODE = Eq(p.diff(t),-2*p+1) pprint(stabODE) # solve an ODE stabSol = dsolve(stabODE,p,hint='best') pprint(stabSol) C1 = Symbol('C1') # substitude C1, so that p(0)=1 stabFun = stabSol.rhs.subs([(C1,Rational(1,2))]) pprint(stabFun) # # plot the solution lam_stab = lambdify(t, stabFun) #set the limits x_min = 0. x_max = 5. y_min = 0.49 y_max = 1. Nt = 100 T = np.linspace(0.,x_max,Nt) P = [lam_stab(Te) for Te in T] fstab = figure(figsize=(12, 12),dpi=120) plot(T, P, lw=2) # generate the direction field X,Y = meshgrid( arange(x_min-0.2,x_max+0.2,0.2),arange(y_min-0.03,y_max+0.03,0.03) ) U = np.copy(X) # adjust the direction field direction to the plot axis U.fill((y_max-y_min)/(x_max-x_min)) V = -2*Y+1 plotDirectionField(X, Y, U, V, "") xlim([x_min,x_max]) ylim([y_min,y_max]) xlabel('t') ylabel('p') show() # TEnd = 5; n = 5; tau = TEnd/n # # We apply the midpoint rule: # pp3 = zeros(n+1) pp3[0] = 1 pp3[1] = lam_stab(tau) lam_stab_rhs = lambdify(x,stabODE.rhs.subs([(p,x)])) for k in range(1,n): # use pnew as the value in the next step pp3[k+1] = pp3[k-1] + 2*tau*lam_stab_rhs(pp3[k]) pprint(pp3) # tp = [(k-1)*tau for k in range(1,n+2)] fstab2 = figure(figsize=(12, 12),dpi=120) plot(T, P, lw=2) plot(tp,pp3,lw=2,c='g') # generate the direction field X,Y = meshgrid( arange(x_min-0.2,x_max+0.2,0.2),arange(y_min-0.03,y_max+0.03,0.03) ) U = np.copy(X) # adjust the direction field direction to the plot axis U.fill((y_max-y_min)/(x_max-x_min)) V = -2*Y+1 plotDirectionField(X, Y, U, V, "") xlim([x_min,x_max]) ylim([y_min,y_max]) xlabel('t') ylabel('p') show() # # Stiff Equation # # define an ODE stiffODE = Eq(p.diff(t),100-100*p) pprint(stiffODE) # solve an ODE stiffSol = dsolve(stiffODE,p,hint='best') pprint(stiffSol) # substitude C1, so that p(0)=2 stiffFun = stiffSol.rhs.subs([(C1,1.0)]) pprint(stiffFun) # # Heun's method fails # #from scipy.integrate import ode lam_stiff_rhs = lambdify(x,stiffODE.rhs.subs([(p,x)])) t0, p0 = 0., 2. TEnd = 1. n = 49 tau = TEnd/n def f(t, p): return lam_stiff_rhs(p) # Heun's method a1 = 0.5 a2 = 0.5 p1 = 1.0 q11 = 1.0 pp4 = zeros(n+1) pp4[0] = p0 for k in range(0,n): tk = t0 + tau*k k1 = f(tk, pp4[k]) k2 = f(tk+p1*tau,pp4[k]+q11*k1*tau) pp4[k+1] = pp4[k] + (a1*k1+a2*k2)*tau ### built-in scipy ode solver # r = ode(f).set_integrator('vodi',method='adams',order=1,nsteps=1) # r.set_initial_value(p0, t0) # while r.successful() and r.t < TEnd: # r.integrate(r.t+tau) # print("%g %g" % (r.t, r.y)) # fstiff = figure(figsize=(12, 12),dpi=120) tp4 = [(k-1)*tau for k in range(1,n+2)] plot(tp4,pp4,lw=2,c='b') #set the limits x_min = 0. x_max = 1. y_min = 0. y_max = 9. # generate the direction field X,Y = meshgrid( arange(x_min-0.05,x_max+0.05,0.05),arange(y_min-0.5,y_max+0.5,0.5) ) U = np.copy(X) # adjust the direction field direction to the plot axis U.fill((y_max-y_min)/(x_max-x_min)) V = 100-100*Y plotDirectionField(X, Y, U, V, "") xlim([x_min,x_max]) ylim([y_min,y_max]) xlabel('t') ylabel('p') show() # # Unless we specify a very small stepsize: # n = 100 tau = TEnd/n pp5 = zeros(n+1) pp5[0] = p0 for k in range(0,n): tk = t0 + tau*k k1 = f(tk, pp5[k]) k2 = f(tk+p1*tau,pp5[k]+q11*k1*tau) pp5[k+1] = pp5[k] + (a1*k1+a2*k2)*tau # tp5 = [(k-1)*tau for k in range(1,n+2)] fstiff2 = figure(figsize=(12, 12),dpi=120) plot(tp4,pp4,lw=2,c='b') plot(tp5,pp5,lw=2,c='g') plotDirectionField(X, Y, U, V, "") xlim([x_min,x_max]) ylim([y_min,y_max]) xlabel('t') ylabel('p') show() # # Maybe it's easier, if we start a little bit later (where the solution is less steep)? # t0 = 0.02 lam_stiff_fun = lambdify(t, stiffFun) newstart = lam_stiff_fun(t0) pprint(newstart) # # Unfortunately not . . . # n = 49 tau = TEnd/n pp6 = zeros(n+1) pp6[0] = newstart for k in range(0,n): tk = t0 + tau*k k1 = f(tk, pp6[k]) k2 = f(tk+p1*tau,pp6[k]+q11*k1*tau) pp6[k+1] = pp6[k] + (a1*k1+a2*k2)*tau # tp6 = [t0+(k-1)*tau for k in range(1,n+2)] fstiff3 = figure(figsize=(12, 12),dpi=120) plot(tp4,pp4,lw=2,c='b') plot(tp5,pp5,lw=2,c='g') plot(tp6,pp6,lw=2,c='r') plotDirectionField(X, Y, U, V, "") xlim([x_min,x_max]) ylim([y_min,y_max]) xlabel('t') ylabel('p') show() show() # # But try t=0.1 as start ... # t0 = 0.1 newstart = lam_stiff_fun(t0) pprint(newstart) # n = 49 tau = TEnd/n pp7 = zeros(n+1) pp7[0] = newstart for k in range(0,n): tk = t0 + tau*k k1 = f(tk, pp7[k]) k2 = f(tk+p1*tau,pp7[k]+q11*k1*tau) pp7[k+1] = pp7[k] + (a1*k1+a2*k2)*tau pprint(pp7)