1:#! /usr/bin/env python 2:# 3:# This is a partial solution to assignment 4 4:# for the course Numerical Methods 2, Acadia University 5:# 6:# Daniel Lemire, Ph.D. 7:# April 2002 8:# GPL v. 3 9:# 10:# Requirements: this file needs Numerical Python to run 11:# 12:# Note: Python and Numerical Python are free software, 13:# see www.python.org. 14:####################################################### 15: 16:from Numeric import * 17:import time 18: 19:########################################################## 20:# section 2.1 21:########################################################## 22:# 23:# 1- First question... solve u'(x)= exp(x)*u(x)+ln(x) numerically. 24:# You can do this with your favorite CAS. I used maxima with the following 25:# command : 26:# ode2(diff(u(x),x)=exp(x)*u(x)+log(x),u(x),x); 27:# (Maxima is free software : http://maxima.sourceforge.net) 28:# to get 29:# u(x)= exp(exp(x))*(integrate(exp(-exp(x))*log(x),x)+C) 30:# 31:# the requirement that u(1)=0 is not so trivial... unless you rewrite 32:# the solution as 33:# 34:# u(x)=exp(exp(x))*integrate(exp(-exp(s))*log(s),s,1,x) 35:# 36:# and then life is good again! 37:# 38:# evaluating u(2) can be done numerically... and we get about 5.97795 39:# 40:# using maxima, the syntax is 41:# 42:# exp(exp(2))*romberg(exp(-exp(x))*log(x),x,1,2),numer; 43:# 44:# We can also get a better idea of what the function look like by 45:# sampling a few points... 46:# 47:# u(1.0)=0 48:# u(1.25)=0.04 49:# u(1.5)=0.24 50:# u(1.75)=1.09 51:# u(1.9) = 2.91 52:# u(1.95) = 4.1 53:# u(2.0)= 6.0 54:# 55:# Basically, we've got exponential growth... It grows very fast... 56:# try figuring out what u(10) is!!! (challenge) 57:# 58:# 2- We've got to do a Runge-Kutta analysis... 59:def f(xx,yy): 60: return exp(xx)*yy+log(xx) 61: 62:def last(x): return x[len(x)-1] 63: 64:def computeRK(f,y,x,dx): 65: k1=dx*f(x,y) 66: k2=dx*f(x+dx/2.0,y+k1/2.0) 67: k3=dx*f(x+dx/2.0,y+k2/2.0) 68: k4=dx*f(x+dx,y+k3) 69: return y+k1/6.0+k2/3.0+k3/3.0+k4/6.0 70: 71:def doRungeKutta(f,y,x,dx): 72: y.append(computeRK(f,last(y),last(x),dx)) 73: x.append(last(x)+dx) 74: 75: 76:def solveByRungeKutta(dx): 77: y=[0.0] 78: x=[1.0] 79: while(last(x) < 2): 80: doRungeKutta(f,y,x,dx) 81: return y 82: 83:def question2_1_2(): 84: print last(solveByRungeKutta(1.0)) 85: print last(solveByRungeKutta(1/2.0)) 86: print last(solveByRungeKutta(1/4.0)) 87: print last(solveByRungeKutta(1/8.0)) 88: print last(solveByRungeKutta(1/16.0)) 89: print last(solveByRungeKutta(1/32.0)) 90: print true 91: ## next few lines are a really, really rough estimation 92: print "order of approximation ", log(abs((last(solveByRungeKutta(1.0)) - true) /(last(solveByRungeKutta(1/2.0))) - true))/log(2.0) 93: print "order of approximation ", log(abs((last(solveByRungeKutta(1/2.0)) - true) /(last(solveByRungeKutta(1/4.0))) - true))/log(2.0) 94: print "order of approximation ", log(abs((last(solveByRungeKutta(1/4.0)) - true) /(last(solveByRungeKutta(1/8.0))) - true))/log(2.0) 95: 96: 97:# 98:# The problem should appear stable. It does appear to converge 99:# fairly fast... the approximation order certainly appears to 100:# be >2 from numerical results. 101:# 102:# 3- The problem with solving the problem over [0,1] is that 103:# f(x,y) blows up at zero. The solution would be to start 104:# at one and stop short of 0, say 0+epsilon. A more involved 105:# solution is to do a change of variable... set w=log(x) 106:# then x= exp(w) and we get a brand new differential 107:# equation where 0 is replaced by -Infinity. Both are 108:# reasonable approaches. 109:# 110:########################################################## 111:# Section 2.2 112:########################################################## 113:# 1- we need to implement Adams-Bashforth 114: 115:def secondLast(x): return x[len(x)-2] 116: 117:def thirdLast(x): return x[len(x)-2] 118: 119:def fourthLast(x): return x[len(x)-3] 120: 121:# Adams-Bashforth 122:def doAdamsBashforth(f,y,x,dx): 123: y3 = last(y) 124: x3 = last(x) 125: y2 = secondLast(y) 126: x2 = secondLast(x) 127: y1 = thirdLast(y) 128: x1 = thirdLast(x) 129: y0 = fourthLast(y) 130: x0 = fourthLast(x) 131: # Something like the next line is really a third order, but was accepted without penalty 132: # because I might not have been clear by what I meant by "order". Conclusions should be the same. 133: #y.append(y2+dx/12.0*(23.0*f(x2,y2)-16.0*f(x1,y1)+5.0*f(x0,y0))) 134: y.append(y3+dx/24.0*(55*f(x3,y3)-59.0*f(x2,y2)+37.0*f(x1,y1)-9*f(x0,y0))) 135: x.append(last(x)+dx) 136: 137:def solveByAdamsBashforth(dx): 138: y=[0.0] 139: x=[1.0] 140: doRungeKutta(f,y,x,dx) 141: doRungeKutta(f,y,x,dx) 142: #print y 143: while(last(x) < 2): 144: doAdamsBashforth(f,y,x,dx) 145: return y 146: 147:def question2_2_1(): 148: print last(solveByAdamsBashforth(1/4.0)) 149: print last(solveByAdamsBashforth(1/8.0)) 150: print last(solveByAdamsBashforth(1/16.0)) 151: print last(solveByAdamsBashforth(1/32.0)) 152: print last(solveByAdamsBashforth(1/64.0)) 153: print last(solveByAdamsBashforth(1/128.0)) 154: print last(solveByAdamsBashforth(1/256.0)) 155: print last(solveByAdamsBashforth(1/512.0)) 156: print last(solveByAdamsBashforth(1/1024.0)) 157: print last(solveByAdamsBashforth(1/2048.0)) 158: print last(solveByAdamsBashforth(1/4096.0)) 159: true = last(solveByAdamsBashforth(1/8192.0)) 160: print true 161: print "order of approximation ", log(abs((last(solveByAdamsBashforth(1.0)) - true) /(last(solveByAdamsBashforth(1/2.0))) - true))/log(2.0) 162: print "order of approximation ", log(abs((last(solveByAdamsBashforth(1/2.0)) - true) /(last(solveByAdamsBashforth(1/4.0))) - true))/log(2.0) 163: print "order of approximation ", log(abs((last(solveByAdamsBashforth(1/4.0)) - true) /(last(solveByAdamsBashforth(1/8.0))) - true))/log(2.0) 164: 165:# This approach is noticeably faster than Runge-Kutta for the same number 166:# of nodes (about 20% with this implementation). It is stable (maybe 167:# a little less stable than RK). With a more efficient implementation, Adams-Bashforth 168:# could be substantially faster. The choice between Runge-Kutta and Adams-Bashforth 169:# could be a difficult one. If you need a lot of nodes, Adams-Bashforth is clearly 170:# faster, so that would be a better choice. However, for this case, Runge-Kutta 171:# appears more accurate so it could be the best choice. 172:# 173:# The lesson here is that Runge-Kutta is not made obselete by multistep approaches. 174:# 175:# 2- we need to implement Adams-Moulton 176: 177:# Adams-Moulton 178:def doAdamsMoulton(f,y,x,dx): 179: # of course, we could optimize this a bit 180: # quite a bit actually! shouldn't call f all that much! 181: y4 = y.pop() 182: x4 = x.pop() 183: y3 = last(y) 184: x3 = last(x) 185: y2 = secondLast(y) 186: x2 = secondLast(x) 187: y1 = thirdLast(y) 188: x1 = thirdLast(x) 189: y0 = fourthLast(y) 190: x0 = fourthLast(x) 191: # Something like the next line is really a fourth order, but was accepted without penalty 192: # because I might not have been clear by what I meant by "order". Conclusions should be the same. 193: # y.append(y2+dx/24.0*(9*f(x3,y3)+19*f(x2,y2)-5*f(x1,y1)+f(x0,y0))) 194: y.append(y3+dx/720.0*(251*f(x4,y4)+646*f(x3,y3)-264*f(x2,y2)+106*f(x1,y1)-19*f(x0,y0))) 195: x.append(last(x)+dx) 196: 197:def solveByAdamsMoulton(dx): 198: y=[0.0] 199: x=[1.0] 200: doRungeKutta(f,y,x,dx) 201: doRungeKutta(f,y,x,dx) 202: #print y 203: while(last(x) < 2): 204: doAdamsBashforth(f,y,x,dx) 205: doAdamsMoulton(f,y,x,dx) 206: return y 207: 208:def question2_2_2(): 209: print last(solveByAdamsMoulton(1/4.0)) 210: print last(solveByAdamsMoulton(1/8.0)) 211: print last(solveByAdamsMoulton(1/16.0)) 212: print last(solveByAdamsMoulton(1/32.0)) 213: print last(solveByAdamsMoulton(1/64.0)) 214: print last(solveByAdamsMoulton(1/128.0)) 215: print last(solveByAdamsMoulton(1/256.0)) 216: print last(solveByAdamsMoulton(1/512.0)) 217: print last(solveByAdamsMoulton(1/1024.0)) 218: print last(solveByAdamsMoulton(1/2048.0)) 219: true = last(solveByAdamsMoulton(1/4096.0)) 220: print true 221: print "order of approximation ", log(abs((last(solveByAdamsMoulton(1.0)) - true) /(last(solveByAdamsMoulton(1/2.0))) - true))/log(2.0) 222: print "order of approximation ", log(abs((last(solveByAdamsMoulton(1/2.0)) - true) /(last(solveByAdamsMoulton(1/4.0))) - true))/log(2.0) 223: print "order of approximation ", log(abs((last(solveByAdamsMoulton(1/4.0)) - true) /(last(solveByAdamsMoulton(1/8.0))) - true))/log(2.0) 224:# 225:# Adams-Moulton is stable and far more accurate than Adams-Bashforth. We might 226:# still prefer Runge-Kutta though in this case! However, if we need a lot 227:# of nodes with ok accuracy, Adams-Moulton might be better. 228:# 229: 230:########################################################## 231:# Section 2.3 232:########################################################## 233:# 234:# The major difference here is that we have 235:# a second order (linear) differential equation. 236:# 237:# We need to reduce this second order equation into 238:# two first order differential equations. 239:# 240:# we set u=y and v = u'=y' 241:# 242:# and we get... 243:# u'= v 244:# v'=exp(x)*v+log(x) 245:# with initial conditions 246:# u(0)=0 247:# and 248:# v(0)= t where t is unknown (at first) 249:# 250:# The problem didn't specify which method we had to use... 251:# a clever student would use Euler!!! 252:# a really smart student would read on the topic and use 253:# an advanced approach... we go with Euler for now... 254: 255:def f1(x,u,v): 256: return v 257: 258:def f2(x,u,v): 259: return exp(x)*v+log(x) 260: 261: 262:def computeEuler1(f1,u,v,x,dx): 263: return u+dx*f1(x,u,v) 264: 265:def computeEuler2(f2,u,v,x,dx): 266: return v+dx*f2(x,u,v) 267: 268:def doSecondOrder(f1,f2,u,v,x,dx): 269: lastu = last(u) 270: lastv = last(v) 271: u.append(computeEuler1(f1,lastu,lastv,last(x),dx)) 272: v.append(computeEuler2(f2,lastu,lastv,last(x),dx)) 273: x.append(last(x)+dx) 274: 275: 276:def shootSecondOrder(t,dx): 277: u=[0.0] 278: v=[t] 279: x=[1.0] 280: while(last(x) < 2): 281: doSecondOrder(f1,f2,u,v,x,dx) 282: return u 283: 284:def solveSecondOrder(dx): 285: u0=shootSecondOrder(0,dx) 286: u1=shootSecondOrder(1,dx) 287: a=last(u0) 288: b=last(u1) 289: if ((a-b) == 0): 290: print "we're dead" 291: else: 292: alpha = (pi - b)/(a-b) 293: solution = alpha*array(u0)+(1-alpha)*array(u1) 294: print "dx = ", dx 295: print "y(0) = ", solution[0] 296: print "y(1.25) = ", solution[len(solution)/4] 297: print "y(1.5) = ", solution[len(solution)/2] 298: print "y(1.75) = ", solution[3*len(solution)/4] 299: print "y(2) = ",last(solution) 300: 301:def question2_3(): 302: solveSecondOrder(1.0) 303: solveSecondOrder(1.0/2.0) 304: solveSecondOrder(1.0/4.0) 305: solveSecondOrder(1.0/8.0) 306: solveSecondOrder(1.0/16.0) 307: solveSecondOrder(1.0/32.0) 308: solveSecondOrder(1.0/64.0) 309: solveSecondOrder(1.0/128.0) 310: solveSecondOrder(1.0/256.0) 311: solveSecondOrder(1.0/512.0) 312: solveSecondOrder(1.0/1024.0) 313: solveSecondOrder(1.0/2048.0) 314: solveSecondOrder(1.0/4096.0) 315: 316:# 317:# You were expected to plot y(1.25), y(1.5) and y(1.75). You should 318:# verify that y(2) is indeed pi all the time...! 319:# 320: 321:# 322:# The following function is just a general call to all of the answers! 323:# 324:if __name__ == '__main__': 325: before = time.time() 326: print "Question 2.1.2" 327: question2_1_2() 328: after = time.time() 329: print "Time used up : ", after-before 330: print "Question 2.2.1" 331: before = time.time() 332: question2_2_1() 333: after = time.time() 334: print "Time used up : ", after-before 335: print "Question 2.2.2" 336: before = time.time() 337: question2_2_2() 338: after = time.time() 339: print "Time used up : ", after-before 340: print "Question 2.3" 341: question2_3() 342: 343: 344: 345: 346: 347: 348: 349: 350: 351: