(*^ ::[ Information = "This is a Mathematica Notebook file. It contains ASCII text, and can be transferred by email, ftp, or other text-file transfer utility. It should be read or edited using a copy of Mathematica or MathReader. If you received this as email, use your mail application or copy/paste to save everything from the line containing (*^ down to the line containing ^*) into a plain text file. On some systems you may have to give the file a name ending with ".ma" to allow Mathematica to recognize it as a Notebook. The line below identifies what version of Mathematica created this file, but it can be opened using any other version as well."; FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2"; MacintoshStandardFontEncoding; 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, "Times"; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, 10, "Times"; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-2, 12, "Courier"; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-2, 12, "Courier"; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R32768, L-2, 12, "Courier"; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-2, 12, "Courier"; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B32768, L-2, 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 = leftheader, inactive, L2, 12, "Times"; fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7, 12, "Times"; fontset = leftfooter, inactive, L2, 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"; paletteColors = 128; currentKernel; ] :[font = title; inactive; preserveAspect; startGroup] Linear Methods of Applied Mathematics An example of an integral equation ;[s] 3:0,1;38,2;39,0;74,-1; 3:1,25,18,Times,1,24,0,0,0;1,24,17,Times,1,23,0,0,0;1,12,8,Times,1,11,0,0,0; :[font = text; inactive; preserveAspect; plain; bold] (c) Copyright 1994,1997,2000 by Evans M. Harrell II and James V. Herod. All rights reserved. :[font = text; inactive; Cclosed; preserveAspect; plain; bold; startGroup] Notes for the instructor. ;[s] 2:0,1;25,0;26,-1; 2:1,13,9,Times,1,12,0,0,0;1,11,8,Times,1,9,0,0,0; :[font = text; inactive; preserveAspect; endGroup] This contains calculations and examples which correlate with chapter 14 of the WWW text by Harrell and Herod. Students can be encouraged to cut and paste from this notebook to do homework. :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Instructions :[font = text; inactive; preserveAspect; endGroup] This notebook uses Mathematica to perform calculations for Harrell and Herod's hypertext book, Linear Methods of Applied Mathematics. The student needs only a basic knowledge of Mathematica to use the notebook, which is designed both to show how to work problems in the text and to provide a template for using Mathematica to work other problems of the student's own design. Calculations will be performed when the reader presses enter in a given calculation cell (bold print). It is best to activate the cells in order, so that Mathematica will be able to call on operators and functions defined in earlier cells. Red color coding is used to warn the reader when a given calculation relies on an earlier one. ;[s] 13:0,0;19,1;30,0;96,1;133,0;181,1;192,0;315,1;326,0;538,1;549,0;626,2;642,0;721,-1; 3:7,13,9,Times,0,12,0,0,0;5,13,9,Times,2,12,0,0,0;1,13,9,Times,0,12,65535,0,0; :[font = section; inactive; Cclosed; preserveAspect; startGroup] Example of an Equation with a Small Integral Operator :[font = text; inactive; preserveAspect] Let's consider the integral equation y(x) = K y + f(x) where K g denotes the integral from 0 to x of Sin[Pi(x-t)] g(t) dt, we take f(x) = x, and the variable x runs from 0 to 1. This kernel is not separable, even though a trigonometric identity can be used to write Sin[Pi(x-t)] as a sum in the form ap(x) bp(t), because the limits of the integral depend on x. (The equation is said to be of Volterra type.) We first ask whether this integral kernel is small in the two senses we have learned about. This means checking whether either of two calculations is less than 1. Notice that since the integral runs only from 0 to x, the integral kernel is K(x,t) = Sin[Pi(x-t)] , 0 < t < x < 1 K(x,t) = 0 , 0 < x < t < 1 The integrals we shall calculate will thus only run from 0 to x in the variable t. To see if the kernel is small in the sense appropriate for a continuous solution y(x), we calculate: ;[s] 9:0,0;317,2;318,0;323,2;324,0;409,1;417,0;901,1;912,0;942,-1; 3:5,13,9,Times,0,12,0,0,0;2,13,9,Times,2,12,0,0,0;2,21,13,Times,64,12,0,0,0; :[font = input; preserveAspect] Integrate[Sin[Pi(x-t)], {t,0,x}] :[font = text; inactive; preserveAspect] This is maximized by 2/š, which is less than 1: :[font = input; preserveAspect] N[2/Pi] :[font = text; inactive; preserveAspect] We are thus guaranteed that the integral equation has a continuous solution, and that we can calculate it by iteration. To see if the kernel is small in the sense appropriate for a square-integrable solution y(x), we calculate: ;[s] 3:0,0;182,1;200,0;230,-1; 2:2,13,9,Times,0,12,0,0,0;1,13,9,Times,2,12,0,0,0; :[font = input; preserveAspect] Integrate[Sin[Pi(x-t)]^2, {t,0,x}] :[font = input; preserveAspect] Integrate[%, {x,0,1}] :[font = input; preserveAspect] N[%] :[font = text; inactive; preserveAspect; endGroup] This is also less than 1, so we are guaranteed that iteration will lead us to the solution in the r.m.s. sense. Normally this is a weaker statement than the previous one, so we are not surprised. Next we calculate the solution :[font = section; inactive; Cclosed; preserveAspect; startGroup] Solving the Equation :[font = input; preserveAspect] Clear[K] :[font = input; preserveAspect] K[f_] := Integrate[Sin[Pi(x-t)] (f /. x -> t), {t,0,x}] :[font = subsubsection; inactive; Cclosed; preserveAspect; startGroup] Explanation of the syntax :[font = text; inactive; preserveAspect; endGroup] The possibly puzzling command embedded in this, (f /. x -> t), has the effect of changing the variable in the function f from x to t, which we need to do to integrate properly. :[font = input; Cclosed; preserveAspect; startGroup] y1[x_] = K[x] ;[s] 3:0,0;9,1;10,0;14,-1; 2:2,12,10,Courier,1,12,0,0,0;1,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] x/Pi - Sin[Pi*x]/Pi^2 ;[o] x Sin[Pi x] -- - --------- Pi 2 Pi :[font = input; Cclosed; preserveAspect; startGroup] y2[x_] = K[y1[x]] ;[s] 5:0,0;9,1;10,0;11,1;13,0;18,-1; 2:3,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] (4*Pi*x + 2*Pi*x*Cos[Pi*x] - Sin[Pi*x])/(4*Pi^3) - (5*Sin[Pi*x])/(4*Pi^3) ;[o] 4 Pi x + 2 Pi x Cos[Pi x] - Sin[Pi x] 5 Sin[Pi x] ------------------------------------- - ----------- 3 3 4 Pi 4 Pi :[font = input; Cclosed; preserveAspect; startGroup] y3[x_] = K[y2[x]] ;[s] 5:0,0;9,1;10,0;11,1;13,0;18,-1; 2:3,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect; endGroup] (-23*Sin[Pi*x])/(16*Pi^4) + (16*Pi*x + 14*Pi*x*Cos[Pi*x] - 7*Sin[Pi*x] + 2*Pi^2*x^2*Sin[Pi*x])/(16*Pi^4) ;[o] -23 Sin[Pi x] ------------- + (16 Pi x + 14 Pi x Cos[Pi x] - 4 16 Pi 2 2 4 7 Sin[Pi x] + 2 Pi x Sin[Pi x]) / (16 Pi ) :[font = input; preserveAspect; endGroup] {Plot[x, {x,0,1}], Plot[y1[x]+x, {x,0,1}], \ Plot[y2[x]+y1[x]+x, {x,0,1}], \ Plot[y3[x]+y2[x]+y1[x]+x, {x,0,1}]} :[font = section; inactive; Cclosed; preserveAspect; startGroup] Error Analysis :[font = text; inactive; preserveAspect] These last two plots look somewhat similar, which is an indication that we are getting close to the exact solution. But how far away are we really? Well, if you stop the series at third order, f + Kf + K2 f + K3 f, then the error is K4 f + K5 f +... We know that the calculations we must check for convergence are actually bounds of the form ||K g|| <= c || g||, where for the infinity norm (for continuous functions), c = 2/š, and for the r.m.s. norm (for square-integrable functions) c = 1/2. (The calculation we did gave 1/4, but that is actually c2). The right scale for this problem is ||x||infinity = 1 or, respectively, ||x||2 = 1/Sqrt[3] ;[s] 15:0,0;214,1;215,0;221,1;222,0;247,1;248,0;254,1;255,0;568,1;569,0;616,2;624,0;654,2;655,0;668,-1; 3:8,13,9,Times,0,12,0,0,0;5,21,13,Times,32,12,0,0,0;2,21,13,Times,64,12,0,0,0; :[font = input; preserveAspect] N[1/Sqrt[3]] :[font = text; inactive; preserveAspect] The norm of Kn f is at most Kn times the norm of f (either norm), so the size of the error is c4 + c5 + ... = c4 (1 + c + c2 + ...) = c4/(1-c). Conclusion: The maximum possible error |x+y1+y2+y3 - yexact(x)| is at most ;[s] 17:0,0;13,1;14,0;30,1;31,0;97,1;98,0;102,1;103,0;113,1;114,0;125,1;126,0;137,1;138,0;219,2;224,0;241,-1; 3:9,13,9,Times,0,12,0,0,0;7,21,13,Times,32,12,0,0,0;1,21,13,Times,64,12,0,0,0; :[font = input; preserveAspect] N[(2/Pi)^4 /(1-2/Pi)] :[font = text; inactive; preserveAspect] The maximum possible error in the rms sense is :[font = input; preserveAspect] (1/2)^4 / (1-1/2) :[font = text; inactive; preserveAspect] Which is fairly small. If we want more accuracy, we need more terms. With 8 terms, :[font = input; preserveAspect] N[(2/Pi)^8 /(1-2/Pi)] :[font = input; preserveAspect] N[(1/2)^7] :[font = text; inactive; preserveAspect] So we would be within 5% uniformly and about 1% in the rms sense. Actually, we are probably doingmuch better, for already at third order, :[font = input; preserveAspect] Plot[y3[x], {x,0,1/2}] Plot[y3[x], {x,1/2,1}] ;[s] 5:0,0;5,1;7,0;28,1;30,0;46,-1; 2:3,12,10,Courier,1,12,0,0,0;2,12,10,Courier,1,12,65535,0,0; :[font = input; Cclosed; preserveAspect; startGroup] y4[x_] = K[y3[x]] y5[x_] = K[y4[x]] y6[x_] = K[y5[x]] y7[x_] = K[y6[x]] y8[x_] = K[y7[x]] ;[s] 11:0,0;11,1;13,0;29,1;31,0;47,1;49,0;65,1;67,0;83,1;85,0;91,-1; 2:6,12,10,Courier,1,12,0,0,0;5,12,10,Courier,1,12,65535,0,0; :[font = output; output; inactive; preserveAspect] (-51*Sin[Pi*x])/(32*Pi^5) + (96*Pi*x + 114*Pi*x*Cos[Pi*x] - 2*Pi^3*x^3*Cos[Pi*x] - 57*Sin[Pi*x] + 24*Pi^2*x^2*Sin[Pi*x])/(96*Pi^5) ;[o] -51 Sin[Pi x] ------------- + (96 Pi x + 114 Pi x Cos[Pi x] - 5 32 Pi 3 3 2 2 2 Pi x Cos[Pi x] - 57 Sin[Pi x] + 24 Pi x Sin[Pi x] 5 ) / (96 Pi ) :[font = output; output; inactive; preserveAspect] (-443*Sin[Pi*x])/(256*Pi^6) + (768*Pi*x + 1122*Pi*x*Cos[Pi*x] - 36*Pi^3*x^3*Cos[Pi*x] - 561*Sin[Pi*x] + 282*Pi^2*x^2*Sin[Pi*x] - 2*Pi^4*x^4*Sin[Pi*x])/(768*Pi^6) ;[o] -443 Sin[Pi x] -------------- + (768 Pi x + 1122 Pi x Cos[Pi x] - 6 256 Pi 3 3 36 Pi x Cos[Pi x] - 561 Sin[Pi x] + 2 2 4 4 6 282 Pi x Sin[Pi x] - 2 Pi x Sin[Pi x]) / (768 Pi ) :[font = output; output; inactive; preserveAspect] (-949*Sin[Pi*x])/(512*Pi^7) + (7680*Pi*x + 13110*Pi*x*Cos[Pi*x] - 570*Pi^3*x^3*Cos[Pi*x] + 2*Pi^5*x^5*Cos[Pi*x] - 6555*Sin[Pi*x] + 3660*Pi^2*x^2*Sin[Pi*x] - 50*Pi^4*x^4*Sin[Pi*x])/(7680*Pi^7) ;[o] -949 Sin[Pi x] -------------- + (7680 Pi x + 13110 Pi x Cos[Pi x] - 7 512 Pi 3 3 5 5 570 Pi x Cos[Pi x] + 2 Pi x Cos[Pi x] - 2 2 6555 Sin[Pi x] + 3660 Pi x Sin[Pi x] - 4 4 7 50 Pi x Sin[Pi x]) / (7680 Pi ) :[font = output; output; inactive; preserveAspect] (-4027*Sin[Pi*x])/(2048*Pi^8) + (92160*Pi*x + 178110*Pi*x*Cos[Pi*x] - 9360*Pi^3*x^3*Cos[Pi*x] + 66*Pi^5*x^5*Cos[Pi*x] - 89055*Sin[Pi*x] + 53370*Pi^2*x^2*Sin[Pi*x] - 1020*Pi^4*x^4*Sin[Pi*x] + 2*Pi^6*x^6*Sin[Pi*x])/ (92160*Pi^8) ;[o] -4027 Sin[Pi x] --------------- + (92160 Pi x + 178110 Pi x Cos[Pi x] - 8 2048 Pi 3 3 5 5 9360 Pi x Cos[Pi x] + 66 Pi x Cos[Pi x] - 2 2 89055 Sin[Pi x] + 53370 Pi x Sin[Pi x] - 4 4 6 6 1020 Pi x Sin[Pi x] + 2 Pi x Sin[Pi x]) / 8 (92160 Pi ) :[font = output; output; inactive; preserveAspect; endGroup] (-8483*Sin[Pi*x])/(4096*Pi^9) + (1290240*Pi*x + 2763810*Pi*x*Cos[Pi*x] - 165690*Pi^3*x^3*Cos[Pi*x] + 1680*Pi^5*x^5*Cos[Pi*x] - 2*Pi^7*x^7*Cos[Pi*x] - 1381905*Sin[Pi*x] + 871920*Pi^2*x^2*Sin[Pi*x] - 20580*Pi^4*x^4*Sin[Pi*x] + 84*Pi^6*x^6*Sin[Pi*x])/(1290240*Pi^9) ;[o] -8483 Sin[Pi x] --------------- + (1290240 Pi x + 2763810 Pi x Cos[Pi x] - 9 4096 Pi 3 3 5 5 165690 Pi x Cos[Pi x] + 1680 Pi x Cos[Pi x] - 7 7 2 Pi x Cos[Pi x] - 1381905 Sin[Pi x] + 2 2 4 4 871920 Pi x Sin[Pi x] - 20580 Pi x Sin[Pi x] + 6 6 9 84 Pi x Sin[Pi x]) / (1290240 Pi ) :[font = text; inactive; preserveAspect] What beautiful functions! You can hardly see the difference between the really accurate solution and the third order solution: :[font = input; preserveAspect] Plot[{x+y1[x]+y2[x]+y3[x]+y4[x]+y5[x]+y6[x] \ + y7[x]+y8[x], \ x+y1[x]+y2[x]+y3[x]}, {x,0,1}] ;[s] 23:0,0;8,1;10,0;14,1;16,0;20,1;22,0;26,1;28,0;32,1;34,0;38,1;40,0;52,1;54,0;58,1;60,0;74,1;76,0;80,1;82,0;86,1;88,0;106,-1; 2:12,12,10,Courier,1,12,0,0,0;11,12,10,Courier,1,12,65535,0,0; :[font = input; preserveAspect; endGroup; endGroup] Plot[{x+y1[x]+y2[x]+y3[x]+y4[x]+y5[x]+y6[x] \ + y7[x]+y8[x], \ x+y1[x]+y2[x]+y3[x]}, {x,.99,1}] ;[s] 23:0,0;8,1;10,0;14,1;16,0;20,1;22,0;26,1;28,0;32,1;34,0;38,1;40,0;52,1;54,0;58,1;60,0;74,1;76,0;80,1;82,0;86,1;88,0;108,-1; 2:12,12,10,Courier,1,12,0,0,0;11,12,10,Courier,1,12,65535,0,0; ^*)