Jeremy Knight Math 639 5/4/11 Algorithm 12.3: Crank-Nicolson
<Text-field style="Heading 1" layout="Heading 1">CrankNicolson</Text-field> To approximate the solution to the parabolic partial differential equation LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzZBLUkmbWZyYWNHRiQ2KC1JI21pR0YkNiVRKCYjOTQ4O3VGJy8lJ2l0YWxpY0dRJXRydWVGJy8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictRiM2JS1GLzYlUSgmIzk0ODt0RidGMkY1RjJGNS8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGQi8lKWJldmVsbGVkR1EmZmFsc2VGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYvNiVRInhGJ0YyRjUtSSNtb0dGJDYtUSIsRicvRjZRJ25vcm1hbEYnLyUmZmVuY2VHRkcvJSpzZXBhcmF0b3JHRjQvJSlzdHJldGNoeUdGRy8lKnN5bW1ldHJpY0dGRy8lKGxhcmdlb3BHRkcvJS5tb3ZhYmxlbGltaXRzR0ZHLyUnYWNjZW50R0ZHLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRi82JVEidEYnRjJGNUZURlQtRlE2LVEoJm1pbnVzO0YnRlRGVi9GWUZHRlpGZm5GaG5Gam5GXG8vRl9vUSwwLjIyMjIyMjJlbUYnL0Zib0ZccC1JJW1zdXBHRiQ2JS1GLzYlUSgmYWxwaGE7RicvRjNGR0ZULUYjNiQtSSNtbkdGJDYkUSIyRidGVEZULyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJy1GLDYoLUYjNiYtRlE2LVEifkYnRlRGVkZqb0ZaRmZuRmhuRmpuRlxvRl5vL0Zib0Zgby1GX3A2JS1GLzYlUScmIzk0ODtGJ0ZkcEZULUYjNiVGZ3BGMkY1RltxLUYvNiVRInVGJ0YyRjVGVC1GIzYlLUZfcDYlLUYvNiVRKCYjOTQ4O3hGJ0YyRjVGW3JGW3FGMkY1Rj1GQEZDRkVGSC1GUTYtUSI9RidGVEZWRmpvRlpGZm5GaG5Gam5GXG8vRl9vUSwwLjI3Nzc3NzhlbUYnL0Zib0Zbcy1GaHA2JEZdcUZURlBGYnFGYnFGYnFGYnFGYnFGXXMtRlE2LVEiPEYnRlRGVkZqb0ZaRmZuRmhuRmpuRlxvRmpyRlxzRk1GX3MtRi82JVEibEYnRjJGNUZQRmJxRmJxRmJxRmJxRmJxRl1zRl9zRmRvRl9zLUYvNiVRIlRGJ0YyRjVGVA==. subject to boundary conditions LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzY1LUkjbWlHRiQ2JVEidUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUkjbW5HRiQ2JFEiMEYnL0YzUSdub3JtYWxGJy1JI21vR0YkNi1RIixGJ0Y+LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkYvJSpzeW1tZXRyaWNHRkYvJShsYXJnZW9wR0ZGLyUubW92YWJsZWxpbWl0c0dGRi8lJ2FjY2VudEdGRi8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVRInRGJ0YvRjJGPkY+LUZBNi1RIj1GJ0Y+RkQvRkhGRkZJRktGTUZPRlEvRlRRLDAuMjc3Nzc3OGVtRicvRldGW29GKy1GNjYkLUYjNiYtRiw2JVEibEYnRi9GMkZARllGPkY+RmZuRjpGQC1GQTYtUSJ+RidGPkZERmluRklGS0ZNRk9GUUZTL0ZXRlVGZG9GZG9GZG9GZG9GOi1GQTYtUSI8RidGPkZERmluRklGS0ZNRk9GUUZqbkZcb0ZZRmhvLUYsNiVRIlRGJ0YvRjJGPg==, and the initial conditions LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYzLUkjbWlHRiQ2JVEidUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInhGJ0YvRjItSSNtb0dGJDYtUSIsRicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUkjbW5HRiQ2JFEiMEYnRkFGQUZBLUY+Ni1RIj1GJ0ZBRkMvRkdGRUZIRkpGTEZORlAvRlNRLDAuMjc3Nzc3OGVtRicvRlZGW28tRiw2JVEiZkYnRi9GMi1GNjYkLUYjNiRGOkZBRkFGPS1GPjYtUSJ+RidGQUZDRmluRkhGSkZMRk5GUEZSL0ZWRlRGZG9GZG9GZG9GZG9GWC1GPjYtUSYmbGVxO0YnRkFGQ0ZpbkZIRkpGTEZORlBGam5GXG9GOkZoby1GLDYlUSJsRidGL0YyRkE=. CrankNicolson:=proc(f,l,T,alpha,m::integer,N::integer) #l = endpoint #T = maximum time #alpha = constant from equation #m = number of x subintervals #k = number of time subintervals local i,j,h,k,lambda,u,w,L,t,x,z: #option trace: w:=Array(1..m): L:=Array(1..m): h:=l/m: k:=T/N: lambda:=alpha^2*k/h^2: w[m]:=0: #set Initial Values for i from 1 to m-1 do w[i]:=f(i*h): end do: #Solve a tridiagonal linear system using Algorithm 6.7 L[1]:=1.0+lambda: u[1]:=-lambda/(2*L[1]): for i from 2 to m-2 do L[i]:=1+lambda+lambda*u[i-1]/2: u[i]:=-lambda/(2*L[i]): end do: L[m-1]:=1+lambda+lambda*u[m-2]/2: for j from 1 to N do #Current Time t:=j*k: z[1]:=((1-lambda)*w[1]+lambda/2*w[2])/L[1]: for i from 2 to m-1 do z[i]:=((1-lambda)*w[i]+ (lambda/2)*(w[i+1]+w[i-1]+z[i-1]))/L[i]: end do: w[m-1]:=z[m-1]: for i from m-2 to 1 by -1 do w[i]:=z[i]-u[i]*w[i+1]: end do: #Note: t=t_j if j=N then print("t",t): for i from 1 to m-1 do x:=i*h: #Note: w[i]=w_i,j print (x,w[i]): end do: end if: end do: end proc:
f:=x->evalf(sin(Pi*x)); l:=1.0; T:=50.0; Zio2I0kieEc2IkYlNiRJKW9wZXJhdG9yR0YlSSZhcnJvd0dGJUYlLUkmZXZhbGZHJSpwcm90ZWN0ZWRHNiMtSSRzaW5HNiRGK0koX3N5c2xpYkdGJTYjKiZJI1BpR0YrIiIiOSRGNEYlRiVGJQ== JCIjNSEiIg== JCIkKyYhIiI= CrankNicolson(f,1.0,.5,1.0,10,50); NiRRInQ2IiQiKysrKytdISM1 NiQkIisrKysrNSEjNSQiK2BMNzBCISM3 NiQkIisrKysrPyEjNSQiK3NeZyVRJSEjNw== NiQkIisrKysrSSEjNSQiKylHIipbLichIzc= NiQkIisrKysrUyEjNSQiKystVyU0KCEjNw== NiQkIisrKysrXSEjNSQiK3dlYGZ1ISM3 NiQkIisrKysrZyEjNSQiKy8tVyU0KCEjNw== NiQkIisrKysrcSEjNSQiKyZIIipbLichIzc= NiQkIisrKysrISkhIzUkIit5XmclUSUhIzc= NiQkIisrKysrISohIzUkIitjTDcwQiEjNw== JSFH