(*^ ::[paletteColors = 128; currentKernel; fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8, 24, "Times"; ; fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6, 18, "Times"; ; fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6, 14, "Times"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20, 18, "Times"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15, 14, "Times"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12, 12, "Times"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "New York"; ; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "New York"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, 12, "Courier"; ; fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, 10, "Geneva"; ; fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = Left Header, inactive, 12, "Times"; ; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; ; fontset = Left Footer, inactive, 12, "Times"; ; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; ; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 12, "Times"; ;] :[font = text; inactive; locked; preserveAspect; center; ] Numerical Methods: Mathematica Notebooks (c) John H. Mathews, 1996 To accompany the text: Numerical Methods: for Mathematics, Science, and Engineering, 2nd Ed, 1992 Prentice Hall, Inc. Englewood Cliffs, New Jersey 07632 U.S.A. Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6 Prentice Hall, International Editions: ISBN 0-13-625047-5 This free software is compliments of the author. E-mail address: in%"mathews@fullerton.edu" http://titan.fullerton.edu/~mathews/ ;[s] 6:0,1;67,0;92,1;157,3;159,1;232,2;486,-1; 4:1,17,12,New York,0,12,0,0,0;3,23,17,New York,0,18,0,0,0;1,19,14,New York,0,14,0,0,0;1,27,18,New York,32,14,0,0,0; :[font = text; inactive; locked; preserveAspect; center; ] CONTENTS ;[s] 2:0,1;8,0;9,-1; 2:1,17,12,New York,0,12,0,0,0;1,23,17,New York,1,18,0,0,0; :[font = text; inactive; locked; preserveAspect; ] Chapter 10. Solution of Partial Differential Equations Algorithm 10.1 Finite-Difference Solution for the Wave Equation Algorithm 10.2 Forward-Difference Method for the Heat Equation Algorithm 10.3 Crank-Nicholson Method for the Heat Equation Algorithm 10.4 Dirichlet Method for Laplace's Equation ;[s] 2:0,1;54,0;305,-1; 2:1,17,12,New York,0,12,0,0,0;1,17,12,New York,1,12,0,0,0; :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 10.1 (Finite-Difference Solution for the Wave Equation). To approximate the solution for the wave equation utt(x,t) = c2uxx(x,t) over R = {(x,t): 0 <= x <= a, 0 <= t <= b} with u(0,t) = 0, u(a,t) = 0 for 0 <= t <= b and u(x,0) = f(x), ut(x,0) = g(x) for 0 <= x <= a. Section 10.1 Hyperbolic Equations Page 507 ;[s] 11:0,1;66,0;125,4;127,0;136,2;137,0;138,3;140,0;291,3;292,0;331,1;377,-1; 5:5,17,12,New York,0,12,0,0,0;2,17,12,New York,1,12,0,0,0;1,20,13,New York,32,9,0,0,0;2,20,13,New York,64,9,0,0,0;1,20,13,New York,64,10,0,0,0; :[font = text; inactive; locked; dontPreserveAspect; ] h is the step in the x-direction k is the step in the y-direction. :[font = text; inactive; locked; dontPreserveAspect; ] We will solve the problem on the grid u[[i,j]] where i = 1,2,...,n is the extent along the x-axis. and j = 1,2,...,m is the extent along the t-axis. The method proceds by simple iteration. It is known that r must satisfy r <= 1 for the method to be stable. :[font = text; inactive; preserveAspect; ] First, execute the following group of cells: :[font = input; dontPreserveAspect; startGroup; ] Off[General::spell1]; Clear[GridU,n,m,u]; GridU[n_,m_] := Module[{}, Clear[u]; u = Table[1,{n},{m}]; For[i=1,i<=n,i++, u[[i,1]] = f[i]; u[[i,2]] = (1-r^2) f[i] + k g[i] + r^2/2 (f[i+1]+f[i-1]); ]; For[j=1,j<=m,j++, Module[{}, u[[1,j]] = 0; u[[n,j]] = 0; ]]; Print["The set up grid for u[[i,j]] has been formed."];]; :[font = input; dontPreserveAspect; endGroup; ] Clear[SolveU,n,m,u]; SolveU[n_,m_] := Module[{i,j}, For[j=3,j<=m,j++, Module[{}, For[i=2,i<=n-1,i++, u[[i,j]] = (2-2 r^2) u[[i,j-1]] + r^2 (u[[i+1,j-1]]+u[[i-1,j-1]]) - u[[i,j-2]]; ]]]; Print["The solution has been computed."];]; On[General::spell1]; :[font = text; inactive; preserveAspect; ] Chapter 10, Example 10.1, Page 504. Consider the wave equation where c^2 = 1. The length of the rod is L = 1. Assume that the initial position is u(x,0) = f(x) = sin(Pi x) + sin(2Pi x). Compare the solution with the exact solution: u(x,t) = sin(Pi x)cos(Pi t) + sin(2Pi x)cos(4Pi t) Solution. We will use c = 1 h = 0.2. k = 0.04 This forces r = 1^2 0.04/(0.2)^2 = 1. ;[s] 3:0,0;1,1;37,0;457,-1; 2:2,16,12,New York,0,12,0,0,0;1,16,12,New York,1,12,0,0,0; :[font = text; inactive; preserveAspect; ] First enter the boundary conditons: :[font = input; dontPreserveAspect; ] Clear[f,g]; f[i_] := Sin[Pi 1/(n-1) (i-1)] + Sin[2 Pi 1/(n-1) (i-1)]//N; g[i_] := 0.0; :[font = text; inactive; dontPreserveAspect; ] Now set up the table of solutions. :[font = input; dontPreserveAspect; ] c = 2.0; h = 0.1; k = 0.05; r = c k/h; n=11; m=21; GridU[n,m] :[font = input; dontPreserveAspect; ] SolveU[n,m] :[font = input; preserveAspect; ] ListPlot3D[u, AxesLabel->{"t(i)","x(i)","u"},ViewPoint->{4,2,3}]; :[font = text; inactive; preserveAspect; ] Compare with the analytic solution. :[font = input; preserveAspect; ] SetOptions[Plot3D, PlotPoints->{n,m}]; Plot3D[Sin[Pi x]Cos[2 Pi t] + Sin[2 Pi x]Cos[4 Pi t], {x,0,1},{t,0,1}]; :[font = text; inactive; preserveAspect; ] To see the numerical values enter the command: :[font = input; preserveAspect; ] Transpose[N[u,3]] :[font = text; inactive; dontPreserveAspect; ] Chapter 10, Example 10.2, Page 506 Consider the wave equation where c^2 = 1. The length of the rod is L = 1. Assume that the initial position is u(x,0) = f(x) = 5x/3 for 0 < x < 3/5 5/2 - 5x/2 for 3/5 < x < 1 Solution. We will use c = 1 h = 0.2. k = 0.04 This forces r = 1^2 0.04/(0.2)^2 = 1. ;[s] 3:0,0;1,1;36,0;403,-1; 2:2,16,12,New York,0,12,0,0,0;1,16,12,New York,1,12,0,0,0; :[font = text; inactive; dontPreserveAspect; ] Boundary conditons. We need 11 function values for f. :[font = input; dontPreserveAspect; ] f[1] = 0.0; f[2] = 0.1; f[3] = 0.2; f[4] = 0.3; f[5] = 0.4; f[6] = 0.5; f[7] = 0.6; f[8] = 0.45; f[9] = 0.3; f[10] = 0.15; f[11] = 0.0; g[i_] := 0.0; :[font = text; inactive; dontPreserveAspect; ] Now set up the table of solutions. :[font = input; dontPreserveAspect; ] c = 2.0; h = 0.1; k = 0.05; r = c k/h; n=11; m=21; GridU[n,m] :[font = input; dontPreserveAspect; ] SolveU[n,m] :[font = input; preserveAspect; ] ListPlot3D[u, AxesLabel->{"t(i)","x(i)","u"},ViewPoint->{4,2,3}]; :[font = text; inactive; preserveAspect; ] To see the numerical values enter the command: :[font = input; preserveAspect; ] Transpose[N[u,3]] :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 10.2 (Forward-Difference Method for the Heat Equation). To approximate the solution for the heat equation ut(x,t) = c2uxx(x,t) over R = {(x,t): 0 <= x <= a , 0 <= t <= b} with u(x,0) = f(x) for 0 <= x <= a and u(0,t) = c1, u(a,t) = c2 for 0 <= t <= b. Section 10.2 Parabolic Equations Page 516 ;[s] 13:0,1;66,0;125,2;126,0;135,3;136,0;137,2;139,0;260,2;261,0;276,2;277,0;300,1;345,-1; 4:6,16,12,New York,0,12,0,0,0;2,16,12,New York,1,12,0,0,0;4,19,13,New York,64,9,0,0,0;1,19,13,New York,32,9,0,0,0; :[font = text; inactive; locked; dontPreserveAspect; ] We will solve the problem on the grid u[[i,j]] where i = 1,2,...,n is the extent along the x-axis. and j = 1,2,...,m is the extent along the t-axis. The function g1[j] furnishes the values u[[1,j]] and the function g2[j] furnishes the values u[[n,j]] and the function f[i] furnishes the initial condition u[[i,1]]. For some problems the corners must be modified. The grid u[[i,j] is formed with the subroutine GridU[n,m]. :[font = text; inactive; preserveAspect; ] First, execute the following group of cells: :[font = input; dontPreserveAspect; startGroup; ] Off[General::spell1]; GridU[n_,m_] := Module[{}, Clear[u]; u = Table[1,{n},{m}]; For[i=1,i<=n,i++, u[[i,1]] = f[i] ]; For[j=1,j<=m,j++, Module[{}, u[[1,j]] = g1[j]; u[[n,j]] = g2[j]; ]]; Print["The grid u[[i,j]] has been formed."];]; :[font = input; dontPreserveAspect; endGroup; ] SolveU[n_,m_] := Module[{i,j}, For[j=2,j<=m,j++, For[i=2,i<=n-1,i++, u[[i,j]]=(1-2r)*u[[i,j-1]] + r*(u[[i-1,j-1]] + u[[i+1,j-1]]) ];]; Print["The solution has been computed."]; ]; On[General::spell1]; :[font = text; inactive; dontPreserveAspect; ] Chapter 10, Example 10.3, Page 511. Consider the heat equation where c^2 = 1. The length of the rod is L = 1. Assume that the ends of the rod are held at the temperature u=0. Assume that the initial temperature distribution is u(x,0) = f(x) = 4x - 4x^2. Apply the forward difference method with r=1 and obtain temperature distributions for t=0.04, 0.08, 0.12,...,0.40. Solution. We will use c = 1 h = 0.2. k = 0.04 This forces r = 1^2 0.04/(0.2)^2 = 1. ;[s] 3:0,0;1,1;37,0;544,-1; 2:2,16,12,New York,0,12,0,0,0;1,16,12,New York,1,12,0,0,0; :[font = text; inactive; dontPreserveAspect; ] Enter the boundary conditons: :[font = input; dontPreserveAspect; ] Clear[g1,g2,f]; g1[j_] := 0.0; g2[j_] := 0.0; f[i_] := 4*(0.2*(i-1))*(1 - 0.2*(i-1))//N; :[font = text; inactive; dontPreserveAspect; ] Now set up the table of solutions. :[font = input; dontPreserveAspect; ] n=6; m=11; GridU[n,m] :[font = text; inactive; dontPreserveAspect; ] Next, solve it. :[font = input; dontPreserveAspect; ] v = u; c = 1; h = 0.2; k = 0.02; r = c^2 k/h^2; SolveU[n,m] :[font = input; preserveAspect; ] ListPlot3D[u, AxesLabel->{"t(i)","x(i)","u "},ViewPoint->{4,2,3}]; TableForm[N[Transpose[u],6]] :[font = text; inactive; preserveAspect; ] Now investigate what happens when the step size is too large. :[font = input; dontPreserveAspect; ] n=6; m=11; GridU[n,m] :[font = text; inactive; dontPreserveAspect; ] Next, solve it. :[font = input; dontPreserveAspect; ] v = u; c = 1; h = 0.2; k = 1/30; r = c^2 k/h^2; SolveU[n,m] :[font = input; preserveAspect; ] ListPlot3D[u, AxesLabel->{"t(i)","x(i)","u "},ViewPoint->{4,2,3}]; TableForm[N[Transpose[u],6]] :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 10.3 (Crank-Nicholson Method for the Heat Equation). To approximate the solution for the heat equation ut(x,t) = c2uxx(x,t) over R = {(x,t): 0 <= x <= a , 0 <= t <= b} with u(x,0) = f(x) for 0 <= x <= a and u(0,t) = c1, u(a,t) = c2 for 0 <= t <= b. Section 10.2 Parabolic Equations Page 517 ;[s] 13:0,1;63,0;122,2;123,0;132,3;133,0;134,2;136,0;262,2;263,0;278,2;279,0;303,1;348,-1; 4:6,16,12,New York,0,12,0,0,0;2,16,12,New York,1,12,0,0,0;4,19,13,New York,64,9,0,0,0;1,19,13,New York,32,9,0,0,0; :[font = text; inactive; locked; dontPreserveAspect; ] We will solve the problem on the grid u[[i,j]] where i = 1,2,...,n is the extent along the x-axis. and j = 1,2,...,m is the extent along the t-axis. The function g1[j] furnishes the values u[[1,j]] and the function g2[j] furnishes the values u[[n,j]] and the function f[i] furnishes the initial condition u[[i,1]]. For some problems the corners must be modified. The grid u[[i,j] is formed with the subroutine GridU[n,m]. :[font = text; inactive; preserveAspect; ] First, execute the following group of cells: :[font = input; dontPreserveAspect; startGroup; ] Off[General::spell1]; Clear[GridU,n,m,u]; GridU[n_,m_] := Module[{i,j}, Clear[u]; u = Table[1,{n},{m}]; For[i=1,i<=n,i++, u[[i,1]] = f[i] ]; For[j=1,j<=m,j++, Module[{}, u[[1,j]] = g1[j]; u[[n,j]] = g2[j]; ]]; Print["The grid u[[i,j]] has been formed."];]; :[font = input; dontPreserveAspect; ] Clear[CoeffMat,n]; CoeffMat[n_] := Module[{}, Clear[r]; Va = Table[-1,{n-1}]; Va[[n-1]]=0; Vc = Table[-1,{n-1}]; Vc[[1]]=0; Vd = Table[2 + 2/r,{n}]; Vd[[1]]=1; Vd[[n]]=1; Print["The coefficient matrix has been formed."];] :[font = input; dontPreserveAspect; ] Clear[TriMat,a,d,c,b,n]; TriMat[a_,d_,c_,b_,n_] := Module[{a0,b0,c0,d0,r,t}, a0=a; b0=b; c0=c; d0=d; Clear[x0]; x0 = Table[0,{n}]; For[r=2,r<=n,r++, Module[{}, t = a0[[r-1]]/d0[[r-1]]; d0[[r]] = d0[[r]] - t c0[[r-1]]; b0[[r]] = b0[[r]] - t b0[[r-1]] ]; ]; x0[[n]] = b0[[n]]/d0[[n]]; For[r=n-1,1<=r,r--, x0[[r]] = (b0[[r]] - c0[[r]] x0[[r+1]])/ d0[[r]] ]; Return[x0] ]; :[font = input; dontPreserveAspect; endGroup; ] Clear[SolveU,n,m]; SolveU[n_,m_] := Module[{i,j}, Clear[b,x]; b = Table[0,{n}]; x = Table[0,{n}]; For[j=2,j<=m,j++, Module[{}, b[[1]] = g1[j]; b[[n]] = g2[j]; For[i=2,i<=n-1,i++, b[[i]] = u[[i-1,j-1]] + (2/r -2)*u[[i,j-1]] + u[[i+1,j-1]] ]; x = TriMat[Va,Vd,Vc,b,n]; For[i=1,i<=n,i++, u[[i,j]] = x[[i]]]; ]; ]; ]; On[General::spell1]; :[font = text; inactive; dontPreserveAspect; ] Chapter 10, Example 10.4, Page 515. Consider the heat equation where c^2 = 1. The length of the rod is L = 1. Assume that the ends of the rod are held at the temperature u=0. Assume that the initial temperature distribution is u(x,0) = f(x) = sin(Pi x). Apply the Crank-Nicholson method with r=1 and obtain temperature distributions for t=0.04, 0.08, 0.12,...,0.40. Compare the solution with the exact solution: u(x,t) = sin(Pi x) exp(-t/Pi^2). Solution. We will use c = 1 h = 0.2. k = 0.04 This forces r = 1^2 0.04/(0.2)^2 = 1. ;[s] 3:0,0;1,1;37,0;622,-1; 2:2,16,12,New York,0,12,0,0,0;1,16,12,New York,1,12,0,0,0; :[font = text; inactive; dontPreserveAspect; ] First set up the boundary conditons: :[font = input; dontPreserveAspect; ] Clear[g1,g2,f]; g1[j_] := 0.0; g2[j_] := 0.0; f[i_] := Sin[Pi 0.1 (i-1)]+Sin[3 Pi 0.1 (i-1)] //N; :[font = text; inactive; dontPreserveAspect; ] Now set up the table of solutions. :[font = input; dontPreserveAspect; ] n=11; m=11; GridU[n,m] :[font = text; inactive; dontPreserveAspect; ] Setting up the tri-diagonal matrx with n rows. Indeed, we could get away with n-2 rows, but the implementation is nice this way. The following matrix will usually use r = 1. :[font = input; dontPreserveAspect; ] CoeffMat[n] :[font = text; inactive; dontPreserveAspect; ] Next, solve it. :[font = input; dontPreserveAspect; ] c = 1; h = 0.1; k = 0.01; r = c^2 k/h^2; SolveU[n,m] :[font = input; preserveAspect; ] ListPlot3D[u, AxesLabel->{"t(i)","x(i)","u "},ViewPoint->{4,2,3}]; :[font = text; inactive; preserveAspect; ] Compare with the analytic solution. :[font = input; preserveAspect; ] SetOptions[Plot3D, PlotPoints->{n,m}]; Plot3D[Sin[Pi x]Exp[-0.1 Pi^2 t]+ Sin[3 Pi x]Exp[-0.1 9 Pi^2 t], {t,0,1},{x,0,1},ViewPoint->{4,2,3}]; :[font = text; inactive; locked; preserveAspect; cellOutline; ] Algorithm 10.4 (Dirichlet Method for Laplace's Equation). To approximate the solution of uxx(x,y) + uyy(x,y) = 0 over R = {(x,y): 0 <= x <= a , 0 <= y <= b} with u(x,0) = f1(x), u(x,b) = f2(x) for 0 <= x <= a and u(0,y) = f3(y), u(a,y) = f4(y) for 0 <= y <= b. It is assumed that Dx = Dy = h and that integers n and m exist so that a = nh and b = mh. Section 10.3 Elliptic Equations Page 531 ;[s] 19:0,1;58,0;98,3;100,0;110,3;112,0;199,3;200,0;215,3;216,0;261,3;262,0;277,3;278,0;324,2;325,0;329,2;330,0;402,1;446,-1; 4:9,16,12,New York,0,12,0,0,0;2,16,12,New York,1,12,0,0,0;2,18,13,Symbol,0,12,0,0,0;6,19,13,New York,64,9,0,0,0; :[font = text; inactive; preserveAspect; ] First, execute the following group of cells: :[font = input; dontPreserveAspect; startGroup; ] Off[General::spell1]; Clear[Dirichlet,q,n,m,w]; Dirichlet[q_List,n_,m_,w_]:= Module[{i,j,k}, Err = 1; k = 0; While[(Err > 0.001) && (k <= 70), Module[{}, Err = 0.0; For[j=2,j<=m-1,j++, Module[{}, For[i=2,i<=n-1,i++, Module[{}, Relax = w/4.0 ( u[[i,j+1]] + u[[i,j-1]] + u[[i+1,j]] + u[[i-1,j]] - 4.0 u[[i,j]]); u[[i,j]] = u[[i,j]] + Relax; Err = Max[Err,Abs[Relax]]; ]]]]; Print["Max grid change = ",Err]; k = k+1]]]; :[font = input; dontPreserveAspect; endGroup; ] Clear[Neumann,qList,n,m,w]; Neumann[q_List,n_,m_,w_]:= Module[{i,j,k}, Err = 1; k = 0; While[(Err > 0.001) && (k <= 70), Module[{}, Err = 0.0; For[j=2,j<=m-1,j++, Module[{}, For[i=2,i<=n-1,i++, Module[{}, Relax = w/4.0 ( u[[i,j+1]] + u[[i,j-1]] + u[[i+1,j]] + u[[i-1,j]] - 4.0 u[[i,j]]); u[[i,j]] = u[[i,j]] + Relax; Err = Max[Err,Abs[Relax]]; ]]]]; For[j=2,j<=m-1,j++, Module[{}, Relax = w/4.0 ( 2 u[[n-1,j]] + u[[n,j-1]] + u[[n,j+1]] - 4.0 u[[n,j]]); u[[n,j]] = u[[n,j]] + Relax; Err = Max[Err,Abs[Relax]]]]; Print["Max grid change = ",Err]; k = k+1]]]; On[General::spell1]; :[font = text; inactive; dontPreserveAspect; ] Chapter 10, Example 10.7, Page 529. Solve Laplace's equation over a 9 by 9 grid with boundary conditions Top: 180 Left: 80 Bottom: 20 Right: 0 ;[s] 3:0,0;1,1;37,0;167,-1; 2:2,16,12,New York,0,12,0,0,0;1,16,12,New York,1,12,0,0,0; :[font = text; inactive; dontPreserveAspect; ] First set up the boundary values. :[font = input; dontPreserveAspect; ] Clear[u]; u = Table[ 70, {9},{9}]; For[i=2,i<9,i++, u[[i,1]]=80] For[i=2,i<9,i++, u[[i,9]]=0] For[j=2,j<9,j++, u[[1,j]]=180] For[j=2,j<9,j++, u[[9,j]]=20] u[[1,1]] = (u[[1,2]]+u[[2,1]])/2; u[[1,9]] = (u[[1,8]]+u[[2,9]])/2; u[[9,1]] = (u[[8,1]]+u[[9,2]])/2; u[[9,9]] = (u[[8,9]]+u[[9,8]])/2; :[font = text; inactive; dontPreserveAspect; ] Next, solve it. :[font = input; dontPreserveAspect; ] w = 4/(2+Sqrt[4-(Cos[Pi/8]+Cos[Pi/8])^2])//N; Dirichlet[u,9,9,w] :[font = input; preserveAspect; ] ListPlot3D[Transpose[u], AxesLabel->{"y(i)","x(i)","u "},ViewPoint->{4,2,3}]; TableForm[N[u,4]] :[font = text; inactive; dontPreserveAspect; ] Chapter 10, Example 10.7, Page 529. Solve Laplace's equation with the Neuman boundary condition over a 9 by 9 grid with boundary conditions Top: 180 Left: 80 Bottom: ux = 0 Right: 0 ;[s] 5:0,0;1,2;36,0;194,1;195,0;222,-1; 3:3,16,12,New York,0,12,0,0,0;1,19,13,New York,64,10,0,0,0;1,16,12,New York,1,12,0,0,0; :[font = text; inactive; dontPreserveAspect; ] First set up the boundary values. :[font = input; dontPreserveAspect; ] Clear[u]; u = Table[ 70, {9},{9}]; For[i=2,i<9,i++, u[[i,1]]=80] For[i=2,i<9,i++, u[[i,9]]=0] For[j=2,j<9,j++, u[[1,j]]=180] For[j=2,j<9,j++, u[[9,j]]=20] u[[1,1]] = (u[[1,2]]+u[[2,1]])/2; u[[1,9]] = (u[[1,8]]+u[[2,9]])/2; u[[9,1]]=80; u[[9,9]]=0; For[j=1,j<=9,j++, u[[9,j]] = u[[9,1]]+ (j-1)/8 (u[[9,9]]-u[[9,1]]); ]; :[font = text; inactive; dontPreserveAspect; ] Next, solve it. :[font = input; dontPreserveAspect; ] w = 4/(2+Sqrt[4-(Cos[Pi/8]+Cos[Pi/8])^2])//N; Neumann[u,9,9,w] :[font = input; preserveAspect; ] ListPlot3D[Transpose[u], AxesLabel->{"y(i)","x(i)","u "},ViewPoint->{4,2,3}]; TableForm[N[u,4]] ^*)