legendre my problem vol 1

Legendre My Problem : Volume I PhL 10.10.09 I came into this document trying to solve a simple set of "integral equations" for f(x): !Syntax Error, Idx f(x) P0(x) = 1 // !Syntax Error, Idx f(x) Pn(x) = δn0 !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4...∞ "Legendre My Problem: solve for f(x) " This seemed like a reasonable and simple problem, but turned out to be a real mess. I thought I was arriving at these equations considering the charged spherical half-shell problem (while reading Stakgold pp 173 Vol II), so I wanted of course to solve these equations to find the charge distribution f(x) in George Green style. I now realize these equations are wrong for the half shell, but I nevertheless got interested in Legendre My Problem in its own right. I first review the subject of complete sets of basis functions for ODE's, and in particular for Legendre polynomials Pn(x), and then I make four "attempts" to solve My Problem. Basically my conclusion is that no solution f(x) exists! What I show, however, is that my various limiting numerical methods to find a solution f(x) all fail. We are seeking a function f(x) whose area is 1, but which gives 0 when integrated against all non-constant polynomials of a complete set defined on a different interval. If we allowed as the area was also 0, we could at least conclude that f(x) = 0 and a solution would exist. Intuitively, as n increases in the second equation, you enforce more and more severe conditions on f(x) which try to drive it to f(x) = 0, but the first condition requires the area to be unity. I suspect the same conclusions obtain if you replace the Pn(x) by an arbitrary set of orthogonal polynomials defined on an interval that differs from (0,1). Knowing the "projections" fn of an assumed-extant f(x) on a set of functions which are not complete on the projecting interval does not guarantee that an f(x) exists which has those projections! In our example here we have projections fn = δn0 and indeed, no f(x) exists which has these projections on (0,1). Overview (3.5 pages 9.23.10) 1 1. Background on the Regular Boundary Value Problem 4 2. Application to the Legendre Polynomials. 6 3. Is the Legendre Transform valid for functions of the form f(x)θ(x) ? [yes] 6 4. My Problem, Plan A (attempt to use the basic Legendre Transform -- fails) 9 5. My Problem, Plan B (trying even and odds in memory of the Laplace Transform -- fails) 11 6. My Problem, Plan C (expand on Pn(-1+2x) which form a complete set on 0,1) 12 7. My Problem, Plan D (numeric subdivision of the interval) 23 _____________________________________________________________________________________ Overview (3.5 pages 9.23.10) In Section 1, I review the general Stakgold theory of symmetric differential (2nd order) operators and how this theory leads to a set of eigenfunctions which are orthogonal and complete. I show that any such eigenfunction set implies a certain transform with projection and expansion. It is a good concise review. In Section 2, I apply this general theory to the case of the Legendre operator, and I obtain for my orthonormal and complete eigenfunctions φn(x) = Pn(x) [ (2n+1)/2]1/2. The transform is this: fn = !Syntax Error, Idx f(x) Pn(x) f(x) = (1/2)Σn(2n+1)fn Pn(x) In Section 3 I ask if the transform shown above is valid for f(x) = F(x)θ(x), that is to say, for f(x) which vanishes to the left of x = 0, and then has an arbitrary shape for x ≥ 0. My motivation is that, for a half sphere, one has charge density σ(z) = F(z)θ(z) and this was a problem I had been pondering. Comment: In my Section 1 analysis, I did not state the conditions for which "the transform" is valid, and I don't know where to find this in Stakgold. I think f(x) has to be an element of L2 for the interval of interest, which here is (-1,1), so the condition on f(x) is that it have a finite norm so !Syntax Error, Idx |f(x)|2 < ∞. This will imply that Σn |fn|2 < ∞ by Parseval, which in turn guarantees that all |fn| < ∞ so all projections certainly exist. In this case, I think we can show that the expansion f(x) converges. Note that this does not require that f(x) or its derivatives be continuous or bounded. Example: f(x) = (x - 1/2)-3 is not square integrable on (-1,1). Example: f(x) = 1/x1/3 is square integrable but is not bounded or continuous at the origin. So in terms of f(x) = F(x)θ(x) being "valid" for our Section 2 transform above, we know the answer is YES as long as F(x) is L2 on the interval (0,1). I then consider f(x) = x θ(x) as an example. The projection is then fn = !Syntax Error, Idx x Pn(x) which turns out to vanish for n = 3,5,7... and has fractional values for n = 0,1,2 and all subsequent evens (I list them off). Using Maple, I show how you reconstruct f(x) = Σn(n+1/2)fn Pn(x) by adding an increasing number of terms. I do two more examples with F(x) = 1-x, which has a discontinuity at x = 0, and then F(x) = x2, and do those same Maple plots for each. [ No Maple file exists, but code segments are short and given below.] In the remaining sections, I attend to the main problem of this doc: solving !Syntax Error, Idx F(x) Pn(x) = δn0 In Section 4 (Plan A) I try to "apply" the previous section. I am trying to solve !Syntax Error, Idx F(x) Pn(x) = δn0, so I rewrite this exactly as !Syntax Error, Idx [ F(x)θ(x) ] Pn(x) = δn0 which I regard as a routine Legendre projection. The corresponding expansion then says F(x)θ(x) = 1. But this is wrong for x < 0 since get 0 = 1. The explanation is that if you fully specify the fn (in this case as δn0), then you cannot at the same time partially specify f(x) ( in this case specifying that f(x) = 0 for x < 0) ! I go on to show that if you select f(x) = (1/2)θ(x), you get some other fn 's and then Maple correctly reconstructs f(x). In Section 5 (Plan B) I set out again to solve !Syntax Error, Idx f(x) Pn(x) = δn0 by decomposing f(x) into a sum of even and odd functions. You can do such a separation, but it leads me nowhere. In Section 6 (Plan C) I try yet again to solve !Syntax Error, Idx f(x) Pn(x) = δn0 by using a set of functions that are orthogonal and complete on the restricted interval (0,1). I use φk(x) = (2k+1)1/2 Pk(-1+2x) as an obvious choice and I write down the projection and expansion for these functions. But then when I consider the equation I am trying to solve, I have to expand the Pn(x) in these φk(x) functions as well as f(x), fine: f(x) = Σm fm φm(x) Pn(x) = Σk Fnk φk(x) Fnk = !Syntax Error, Idx φk(x) Pn(x) = known numbers Inserting these expansions into my "equation to be solved" I get the vector equation F f = e0 and my goal is then to solve this for f = F-1e0, then I can construct f(x) = Σm fm φm(x) and I will have the f(x) that solves my problem! To make things more Maple friendly, I shift indices so that Gnm = Fn-1,m-1 and gm = fm-1 so indices on G and g begin with 1 instead of 0. I still have Gg = e1 so g = G-1e0 I then use Maple to compute the matrix G for n = 1 to 6 (N=6), and then the matrix G-1. Then g = G-1e0 is just the first column of G-1. I get these numbers for gm, construct f6(x) = Σm=16 gm φm(x) and plot it in Maple. Indeed this is in fact a solution of my problem !Syntax Error, Idx f(x) Pn(x) = δn0 for n = 0 to 5. The f(x) has a sine shape with a large amplitude on the order of 106. But for n > 5, the integrals do not give δn0= 0. So my idea is then to keep increasing N to ∞, and then I will have my desired final solution! The big question is: what happens as you try to do this? A sort of "experimental approach" to a solution. So I start with N = 6, then I do N = 10, then N = 14. The amplitude of the solution fN(x) gets larger and larger, not a good sign. Maple could not do N = 20, too slow, but I found a way around this and then got to N = 40. Here is a plot of our solution f(x) for N = 40, I show blowups of the two ends just to show that nothing divergent is occurring at either end. A careful count shows 39 zero crossings, as expected. The sine amplitude is now on the order of 1029 !! I could not get it to do N = 100. My conclusion is that fN(x) is not converging to anything, so my method is failing. In my "Conclusion" section, I reexamine my matrix solution fN = FN-1e0 . I have stable values at the start of the sequence fm defined as follows f0 = limN→∞ f(N)0 = 1 f1 = limN→∞ f(N)1 = f2 = limN→∞ f(N)2 = 3 ..... but I conclude that limm→∞ fm = ∞. When we try to expand f(x) = Σm=0∞fm φm(x), we find that the coefficients fm keep getting larger and larger as we go deeper into the series. I am sure this means that the series itself does not converge (perhaps use a comparison test). If our expansion for f(x) does not converge, we are certainly not going to find f(x) by this method. I then show why it is that the fm increase so dramatically (reason 2 is better than reason 1). In any event, in my Final Conclusion, I conclude that fN(x) does not converge to any finite function, though you might vaguely argue that it approaches 0 in a distribution sense. The requirement that the integral of f(x) be 1 is lost in the huge amplitudes involved. In Section 7 (Plan D) I arrive at a similar matrix equation to that obtained above, f = F-1 e1 , but by a different path. I put the integral equation !Syntax Error, Idx f(x) Pn(x) = δn0 on a lattice with N equal intervals, and then fi = f(xi) and Fni = Pn-1(xi)/N, so the meaning of both f and F is different from Section 6 above. For various values of N I have Maple plot the resulting fi points as the first column of F-1. I find that as N increases, the plot of fi has this shape, shown here for N = 40: The pulse moves toward the center of the (0,1) range and gets larger and frequency increases as well. Again, we get the feeling that our solution f(x) is basically f(x) = 0 as a distribution in the limit N → ∞, and again the requirement that the integral of the fi be 1 is completely lost in the huge amplitudes involved. So once again, we conclude that f(x) cannot be found by this method and probably does not exist. __________________________________________________________________________________ 1. Background on the Regular Boundary Value Problem In Stakgold Volume I Section 4.2 we studied the properties of self-adjoint second-order differential operators L. These have a general form L = -∂x(p(x)∂x) + q(x) where functions p and q are real and continuous, and where p is positive in our finite interval (a,b). Unmixed boundary conditions are required at the endpoints to achieve full self-adjointness and then L is "symmetric". We considered the eigenvalue problem Lφλ = λs(x)φλ where we introduced a third function s (called a weighting function) which is also real, continuous and positive. We could have included s(x) in our definition of L, redefining L' = L/s but we did not do this and followed tradition. Thus, it was convenient to define Lλ = L - λs(x) and then our eigenfunctions are the homo solutions of Lλφλ= 0. The ODE system defined in this manner is known as the "regular boundary value problem". One can show that the eigenfunctions φλ of different eigenvalues are orthogonal with weight s(x), and the eigenvalue spectrum is discrete. We then wrote down a Green's Equation of the form Lλg(x|ξ;λ) = δ(x-ξ). Instead of requiring g = 0 at both endpoints (as we do in Volume II for PDE's), we required that Ba,b(g) = 0 (our unmixed BC's) at the two endpoints. We were always able to solve the Green's Equation for the Green's Function, and here was that solution: g(x|ξ; λ) =[1/C(λ)] w(x<, λ) z(x>,λ) C(λ) = p(x) W [ w(x,λ), z(x,λ) ; x ] where here, w(x;λ) satisfies the left end BC and z(x;λ) satisfies the right end BC, and W is the Wronskian. We learned that C(λ) has zeroes at the eigenvalues λn and thus g(x|ξ; λ) has poles at the λn. When we recast our ODE into a symmetric-kernel integral equation, we found that the functions ψn(x) = φn(x) were eigenfunctions of the corresponding integral operator and formed a complete set (with no weighting function) and therefore the φn(x) form a complete set with weight s(x). The complete set idea could be expressed in this manner: Σn φn(x) n(ξ) = δ(x-ξ)/s(x) and !Syntax Error, I dx s(x) φn(x) m(x) = δn,m completeness = closure orthonormality These facts just stated always imply a specific projection/recovery transform as follows: fn = !Syntax Error, Idx s(x) f(x) n(x) f(x) = Σn fn φn(x) That the transform is valid can be shown as follows. First, f(x) = Σn [fn] φn(x) = Σn [!Syntax Error, Idx' s(x')f(x') n(x')] φn(x) = !Syntax Error, Idx' s(x') f(x') {(Σn n(x') φn(x)} = !Syntax Error, Idx' s(x') f(x') {δ(x-x')/s(x')} = !Syntax Error, Idx' f(x') δ(x-x') = f(x) and then, in the other "direction", fn = !Syntax Error, Idx s(x) [f(x)] n(x) = !Syntax Error, Idx s(x) [Σm fm φm(x)] n(x) = Σm fm {!Syntax Error, Idx s(x) φm(x) n(x) } = Σm fm δn,m = fn In words: (1) any reasonable function f(x) on (a,b) can be expanded on the φn(x) with coefficients fn. (2) any set of coefficients fn is associated with some function f(x) on (a,b) according to f(x) = Σn fn φn(x). One could add an arbitrary constant k or kn and then the transform would be fn = k !Syntax Error, Idx s(x) f(x) n(x) f(x) = (1/k) Σn fn φn(x) fn = kn !Syntax Error, Idx s(x) f(x) n(x) f(x) = Σn (1/kn)fn φn(x) 2. Application to the Legendre Polynomials. The ODE is this from Schaum p 146 (add overall minus sign) [ form: Lλ = -∂x(p(x)∂x) + q(x) - λs(x) ] Lλ = –(1-x2)∂x2 + 2x∂x – n(n+1) s(x) = 1 λ = n(n+1) = – ∂x[ (1-x2)∂x] – n(n+1) p(x) = (1-x2) q(x) = 0 Notice that p(x) is non-negative in (-1,1), the interval for this problem. The solutions to this ODE are the eigenfunctions of Lφλ = λφλ and are the Pn(x). Warning: these are orthogonal but not orthonormal. They are normalized so that Pn(1) = 1 instead of as orthonormals. It is easy to form orthonormal versions: φn(x) = Pn(x) [ (2n+1)/2]1/2 Then !Syntax Error, Idx φn(x) φn(x) = [ (2n+1)/2] !Syntax Error, Idx Pn(x) Pn(x) = 1 // Schaum p 147 Our projection/recovery transform can be written fn = kn !Syntax Error, Idx f(x) n(x) f(x) = Σn (1/kn)fn φn(x) fn = kn !Syntax Error, Idx f(x) Pn(x) [ (2n+1)/2]1/2 f(x) = Σn (1/kn)fn Pn(x) [ (2n+1)/2]1/2 If we then select kn = [ (2n+1)/2]-1/2 we get this transform fn = !Syntax Error, Idx f(x) Pn(x) f(x) = (1/2)Σn(2n+1)fn Pn(x) 3. Is the Legendre Transform valid for functions of the form f(x)θ(x) ? [yes] Suppose we restrict our interest to functions f(x) which vanish on (-1,0). Is our Legendre Transform still valid? Let's try a simple example. f(x) = 0 on the left, and f(x) = x on the right. This is a nice example because it is continuous at x = 0. Then fn = !Syntax Error, Idx x Pn(x) How do we compute this integral? This is a bigger problem than I thought, analytically speaking, and that I think lies at the core of my confusion which led me to ponder this problem. The best closed form I can get follows this method: (2n+1) xPn = (n+1)Pn+1 + nPn-1 // Schaum p 147 top line (2n+1) !Syntax Error, Idx xPn = (n+1) !Syntax Error, Idx Pn+1 + n!Syntax Error, Idx Pn-1 = {(n+1) [ (Pn+2 - Pn)/(2n+3) ] + n[ (Pn - Pn-2)/(2n-1) ]} |10 // Schaum p 147 25.34 = - {(n+1) [ (Pn+2(0) - Pn(0))/(2n+3) ] + n[ (Pn(0) - Pn-2(0))/(2n-1) ]} since at x = 1 both terms are 0 (as long as n > 1). Then we find that !Syntax Error, Idx xPn = - {(n+1) [ (Pn+2(0) - Pn(0))/(2n+3) ] + n[ (Pn(0) - Pn-2(0))/(2n-1) ]}/(2n+1) n>1 = - (n+1) /(2n+1) Pn+2(0) + (n+1)/[ (2n+1) (2n+3) ] Pn(0) - (n+1)n/(2n+1) Pn(0) + (n+1) n [(2n-1) (2n+1)] Pn-2(0) For n = 0 we get (1/2) and for n = 1 we get 1/3. Then in Schaum p 147 5.32 we find that Pn(0) = (-1)n/2 (n-1)!! / n!! n = even Pn(0) = 0 n = odd Thus, we learn that for n = 3,5,7...9 fn ≡ !Syntax Error, Idx xPn = 0. And for n even we get a messy expression involving double factorials, I won't bother to write it, but at least it gives us an answer without any sums. Here are some Maple computations. For any specific Pn, Maple of course know how to integrate powers of x to get the answers: fn = !Syntax Error, Idx x Pn(x) f0= 1/2 f1= 1/3 f2= 1/8 f4= -1/48 f6 = 1/128 f8 = -1/256 f10= 7/3072 Now, just using these values, let's try to reconstruct our starting function f(x) = (1/2)Σn(2n+1)fn Pn(x) This is where Maple is going to sparkle! If we keep only f0, f1 and f2 we get this situation: which is not too bad really. Now add in f4 and then f6 and so on to get this pix: 4 6 8 20 The message for me here is that this is essentially a numerical problem! Now suppose we pick the much less pleasant function f(x) = (1-x) instead of x for our positive-only test function. Then here is what we get going out through n = 40 So I think we have an answer to our question. And one other example f(x) = x2 I of course knew the thing was "still valid" but it is nice to see how this works numerically. 4. My Problem, Plan A (attempt to use the basic Legendre Transform -- fails) Solve these "integral equations" for f(x) ( note improper range 0,1 ): !Syntax Error, Idx f(x) P0(x) = 1 !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4... Let's take as our full range function f(x)θ(x). Then suppose we try to do a Legendre expansion as above: f0 = !Syntax Error, Idx f(x)θ(x) P0(x) = !Syntax Error, Idx f(x) P0(x) = !Syntax Error, Idx f(x) = 1 for n = 0 fn = !Syntax Error, Idx f(x)θ(x) Pn(x) = !Syntax Error, Idx f(x) Pn(x) = 0 for n = 1,2,3.... If we insert these projections into our expansion formula, we get f(x)θ(x) = (1/2)Σn=0∞(2n+1)fn Pn(x) = (1/2)f0 = 1/2 but the RHS (namely, 1/2) does not agree with the LHS when x < 0, so something is inconsistent. Also, if we compute the projections of (1/2) θ(x), we do not get 0. For example, !Syntax Error, Idx (1/2) P3(x) = - 1/16. And here is an extra confusion. Suppose we start with f(x) = (1/2)θ(x) and compute its projections, then here are the first 8 fn we obtain: The odd fn are in fact 0, but the evens are not. If we now reconstruct, using 40 terms, we get And now we do in fact reproduce the correct f(x) = (1/2) θ(x). So how do we resolve this mystery? If someone hands me the set of Legendre projections fn = δn,0, there is only one f(x) which generates those projections, and that is f(x) = (1/2) on the full interval (-1,1) It is completely unreasonable to expect that this set of fn would be associated with a function of the form f(x)θ(x). In fact, it is not! If you start with some arbitrary set of fn, you cannot expect to impose conditions on the associated f(x). Here I started with fn = δn,0 so I cannot expect a resulting function of the form f(x)θ(x). So let's review what we learn by applying the Legendre Transform to this problem: (notice full range integrals) f0 = !Syntax Error, Idx f(x) P0(x) = 1 fn = !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4... Solution is f(x) = 1/2. But this is not the problem I am trying to solve, which has 0 to 1 integrals. This shows that you don't learn much by using a set of orthogonal polynomials on the wrong interval! 5. My Problem, Plan B (trying even and odds in memory of the Laplace Transform -- fails) We want to solve these "integral equations" for f(x): !Syntax Error, Idx f(x) P0(x) = 1 !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4... Is there a function f(x) which exactly satisfies this set of equations? Suppose we write f(x) in this way f(x) = fe(x) + fo(x) f(-x) = fe(-x) + fo(-x) = fe(x) – fo(x) which we can solve to get fe(x) = [f(x) + f(-x) ] / 2 fo(x) = [f(x) – f(-x) ] / 2 Then maybe you could say !Syntax Error, Idx f(x) Pn(x) = !Syntax Error, Idx [fe(x) + fo(x)] Pn(x) Now then !Syntax Error, Idx [fe(x)] Pn(x) = !Syntax Error, Idx [fe(x) (-1)n] Pn(x) !Syntax Error, Idx fe(x) Pn(x) = (-1)n!Syntax Error, Idx fe(x)Pn(x) Then !Syntax Error, I dx fe(x) Pn(x) = !Syntax Error, Idx [fe(x)] Pn(x) + (-1)n!Syntax Error, Idx fe(x)Pn(x) = [ 1 + (-1)n] !Syntax Error, Idx[fe(x)] Pn(x) => !Syntax Error, I dx fe(x) Pn(x) = 0 n = odd !Syntax Error, I dx fe(x) Pn(x) = 2 !Syntax Error, Idx[fe(x)] Pn(x) n = even So suppose we assume that our solution f(x) is entirely even, f(x) = fe(x). If that were the case we have => !Syntax Error, I dx f(x) Pn(x) = 0 n = odd !Syntax Error, I dx f(x) Pn(x) = 2 !Syntax Error, Idx[f(x)] Pn(x) = 0 n = 2,4,6,8 !Syntax Error, I dx f(x) P0(x) = 2 !Syntax Error, Idx[f(x)] P0(x) = 1 n = 0 The above three lines say that ALL bona fide Legendre projections are 0 except for n=0, and we are left with f(x) = constant for our answer. But we know this f(x) does not have all zero projections for n>0, so our assumption that f(x) is entirely even is no good. I think if we assumed f(x) was entirely odd, we end up with a similar contradiction. Therefore, if there is a solution f(x), it must be neither even nor odd but must have both "components". 6. My Problem, Plan C (expand on Pn(-1+2x) which form a complete set on 0,1) Solve these "integral equations" for f(x): !Syntax Error, Idx f(x) P0(x) = 1 !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4... Is there a function f(x) which exactly satisfies this set of equations? We know that on the interval (0,1), the Pn(x) are all independent functions since we cannot construct Pn from any lower Pn. The same could be said about the set of powers xn . Our big problem is a lack of orthogonality on this interval. So, we can map 0,1 to -1,1 using the mapping x' = -1+2x. We know that Pn(x) are complete on (-1,1). Then perhaps we can say that Pn(-1+2x) are complete on 0,1. We know that [ (2n+1)/2] !Syntax Error, Idx Pn(x) Pm(x) = δn,m // Schaum p 147 Let x = -1+2x' so that dx = 2dx' so we have (2n+1)!Syntax Error, Idx' Pn(-1+2x') Pm(-1+2x') = δn,m (2n+1)!Syntax Error, Idx Pn(-1+2x) Pm(-1+2x) = δn,m We then have this set of basis functions φn(x) = (2n+1)1/2 Pn(-1+2x) !Syntax Error, Idx φn(x)φm(x) = δn,m Verify: !Syntax Error, Idx φn(x)φm(x) = !Syntax Error, Idx { (2n+1)1/2 Pn(-1+2x)} { (2m+1)1/2 Pm(-1+2x)} = (2n+1)1/2(2m+1)1/2!Syntax Error, Idx Pn(-1+2x) Pm(-1+2x) = (2n+1)1/2(2m+1)1/2δnm/ (2n+1) = δnm OK We think that the projection/transform going with this would be fn = !Syntax Error, Idx φn(x)f(x) f(x) = Σn=0∞fnφn(x) I think the φn(x) form a complete set of basis functions for the interval (0,1). Now go back to !Syntax Error, Idx f(x) P0(x) = 1 !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4... Let's expand (these are new and different coefficients, but I still call them fn) f(x) = Σm fm φm(x) Pn(x) = Σk fk,n φk(x) fk,n = !Syntax Error, Idx φk(x) Pn(x) = known numbers Then we have !Syntax Error, Idx φm(x) Σk fk,n φk(x) = 0 ie !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4... Σm fm Σk fk,n !Syntax Error, Idx φm(x)φk(x) = 0 Σm fm Σk fk,n δm,k = 0 Σm fm fm,n = 0 n = 1,2,3....N-1 Σm fm fm,0 = 1 n = 0 Think of the numbers fm,n as elements of a square matrix of size N, call this thing Fnm . Then Σm=0N-1 Fnm fm = δn,0 n = 0,....N-1 F f = e0 That is, we now assume we have only N terms in our m sum, so that we get square matrices. Assuming det(F) ≠ 0, we can solve for f at least for any finite N according to f = F-1e0 . Fnm = fm,n = !Syntax Error, Idx φm(x) Pn(x) = (2m+1)1/2!Syntax Error, Idx Pm(-1+2x)Pn(x) I am going to use Maple soon, and it does not like matrices indices 0..N-1, it wants indexing to start with 1. So let's define Gnm = Fn-1,m-1 gm = fm-1 G11 = F00 Now consider: Σm=0N-1 Fnm fm = δn,0 Let m = m'-1 and let n = n' - 1 so this says m' = m+1 and n' = n+1 Σm'=1N Fn'-1,m'-1 fm'-1 = δn'-1,0 = δn',1 Σm=1N Fn-1,m-1 fm-1 = δn,1 Σm=1N Gnm gm = δn,1 Gg = e1 g = G-1e1 So we then have Gnm = Fn-1,m-1 = !Syntax Error, Idx φm-1(x) Pn-1(x) = (2m-1)1/2!Syntax Error, Idx Pm-1(-1+2x)Pn-1(x) Now we can construct our matrix G like so, where I pick N = 6 just as an example: [ NOTE: I wrongly wrote 1-2x as P arg instead of -1+2x -- now fixed below! ] If we apply this e1 , we just pick off the first column, so gm = first column Now recall that f(x) = Σm=0N-1 fm φm(x) = Σm=1N fm-1 φm-1(x) = Σm=1N gm φm-1(x) So here I construct the function and check our orthogonality: So here is a function f(x) that exactly meets all our requirements at least for n = 0 through 5. For larger n, we do not get 0 projection for the N=6 calculation, BUT we can increase N to any size we like. For example, if I increase to N = 10, I get these results: (notice how the vertical scale is increasing!) The f(x) is very different now, and it now has the correct projections for the n = 0 through 9. Then as before we have problems. Here is the curve for N = 14 (now vertical extent is in the billions). I have this funny feeling that if we keep increasing N, our f(x) will oscillate faster and faster, and there is no clean limit! The number of zero crossings seems to be N-1. So as N → ∞, we will have an infinite number of zero crossings in the interval (0,1). Also, the amplitude of the waves keeps increasing. Notice, however, that the elements of the solution fm are "stable" as we increase N (compare the above g array to a previous one). This is a characteristic of using an orthonormal set. But the elements of fm seem to be getting larger and larger! Numeric Calculations. Using the method shown above, I was unable to do N = 20 because it ran a full half hour and the inverse matrix did not appear. I then learned that we can do more "numeric" calculations to achieve larger N. The trick is simply to fill the initial matrix with "floating point numbers" like this: G := evalf(matrix(40,40,GG)): Then the remaining code stays the same and things speed up. Here is the plot for N = 40 Notice the numbers are truly huge now, on the order of 1029 . I show blowups of the two ends just to show that nothing divergent is occurring at either end. A careful count shows 39 zero crossings, as expected. The central region is stabilizing into a flat sine pattern. If you now compute the fn = !Syntax Error, Idx f(x) Pn(x) = 0, you no longer get 1,0,0,0...large. Instead, you get small numbers in place of the first 39 0's, and then you get large numbers after that. I did this all with Digits := 20 just for fun. I tried to bump up to N = 100 starting at 1:34. I expect this to increase the time needed to compute the initial matrix G by 2.52 = 6.25 so expecting maybe 5 minutes to do this. But it has run 21 minutes and seems stuck. Memory pins at 116M and the time counter just upticks. The Task Manager does not show a lot of page fault activity. So I am mystified why it is taking so long to build this 100 x 100 matrix. I got it to do 50x50, looks very similar to the above. Conclusion: If we consider these "integral equations", !Syntax Error, Idx fN(x) P0(x) = 1 n = 0 !Syntax Error, Idx fN(x) Pn(x) = 0 n = 1,2,3,4...N there is an exact (and plottable) solution fN(x) for any finite value of N. However, the sequence fN(x) does not seem to converge to anything except possibly some strange distribution that has an infinite number of crossings and has infinite amplitude across the interval! Recall this distributional fact from distributions.doc, 6. limN→∞ sin(Nx) = 0 // I p 44 B If we differentiate this wrt x, we get limn→∞ Ncos(Nx) = 0 This last sequence FN(x) = Ncos(Nx) has properties similar to our fN(x): as N increases, you get more crossings without limit, and the amplitude increases without limit. So one idea is that maybe the distributional limit of our fN(x) is just 0. The only problem is the need for !Syntax Error, Idx f(x) = 1. BUT, since the fn get very large, the function f(x) has a huge "scale" (1029 above!) so the most miniscule alignment of f(x) say near one of the endpoints could easily achieve this area integral of 1 even though f(x) is approaching 0 as a distribution. Yes, this is a fudgy thing to say, but I think it contains the truth of what is going on here. The sequence fm diverges. Let's go back before we got in the above NxN matrix activity. We had g0 ≡ !Syntax Error, Idx f(x) P0(x) = 1 gn ≡ !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4... φn(x) = (2n+1)1/2 Pn(-1+2x) !Syntax Error, Idx φn(x) φm(x) = δn,m fn = !Syntax Error, Idx φn(x)f(x) f(x) = Σn=0∞fnφn(x) Let's expand only f(x) and install into the gn expression n = 1,2,3... gn ≡ 0 = !Syntax Error, Idx Σm=0∞fmφm(x) Pn(x) = Σm=0∞fm {!Syntax Error, Idx φm(x) Pn(x) } = Σm=0∞ Fnm fm OK, this just leads again to our result Σm=0∞ Fnm fm = δn,0 F f = e0 What does this mean as an infinite dimensional matrix? We certainly studied this subject a lot in Stak I. We have merely recast our integral equation into an infinite dimensional matrix problem using our transform. We can then at least use finite matrix techniques and try to find a limit. We have F(N)f(N) = e0(N) f(N) = (F(N))-1 e0(N) and we want to claim that our solution coefficient vector is fm = limN→∞ f(N)m . I think I know that this is true: f0 = limN→∞ f(N)0 = 1 f1 = limN→∞ f(N)1 = f2 = limN→∞ f(N)2 = 3 ..... So we have a sequence here fm . I am sure that limm→∞ fm = ∞. The accumulation point is at ∞. The sequence fm does not converge. When we try to expand f(x) = Σm=0∞fmφm(x), we find that the coefficients fm keep getting larger and larger as we go deeper into the series. I am sure this means that the series itself does not converge (perhaps use a comparison test). If our expansion for f(x) does not converge, we are certainly not going to find f(x) by this method. Is this problem due to my choice of φn functions, or would it be true for any set of functions orthogonal on (0,1), such as sines and cosines ? In any event, we still might claim that f(x) → 0 as a distribution. Question: Why are the fm increasing in size so dramatically. Answer 1 (determinant is large, but not a very good answer). If we go back to our matrix G, we see that the entries get smaller as you move down and to the right. The reason is that in the direction, you are computing the overlap of two highly oscillatory functions and larger n and m means less overlap. Secondly, we get a triangular matrix simply because Pn when expressed in terms of a sum of φm cannot have φm with a power more than xn which means the overlap is 0 when m > n, as we see below. Because the rightmost non-zero number in each row tends to be small, and because for a triangular matrix the determinant is just the product of these small diagonal elements, we get a very small determinant, it being the product of the smallest number on each row. In our example here, det = 10-18 . So that at least provides the potential for an element of the inverse matrix to be large! Now here is the inverse matrix for the above matrix, The first column are the fn in our expansion f(x) = Σn=0∞fnφn(x). Answer 2: (better) Consider this integral and what it is equal to, !Syntax Error, Idx [ f0 + f1φ1(x) + f2φ2(x) + ... + fn φn(x) + ... + f∞ φ∞(x) ] Pn(x) = !Syntax Error, Idx [ f0 + f1φ1(x) + f2φ2(x) + ... + fn φn(x) ] Pn(x) because the higher φm(x) terms m > n have no overlap with Pn(x) (as the triangularity of G verifies). Assume that we have found a set of fi (i = 0..n) such that this last integral is 0 and in fact all "lower" integrals are also 0 except the base integral is 1. As we now add another term fn+1 φn+1(x) to the series, none of the integrals mentioned so far is altered regardless of the value of fn+1! So really our only requirement of this new series term is that we need this to be true: !Syntax Error, Idx [ f0 + f1φ1(x) + f2φ2(x) + ... + fn φn(x) + fn+1 φn+1(x)] Pn+1(x) = 0 which we can write as !Syntax Error, Idx [ f0 + f1φ1(x) + f2φ2(x) + ... + fn φn(x)] Pn+1(x) = – fn+1 !Syntax Error, Idx φn+1(x) Pn+1(x) Now the LHS is presumably a fairly large number. The reason for that is that the scale of f(x) in [..] is itself quite large, since fn is large, and when we move from Pn(x) to Pn+1(x), we cause a misalignment of the two functions ( which were in a sense aligned before so as to give 0 integral) , and so since fn is large, we expect a moderately large number. In our example with n = 9, this number is 1452.3. At the same time, we expect !Syntax Error, Idx φn+1(x) Pn+1(x) to be a small number ( in our case .0002) since it is the overlap of two fast rolling functions of unit amplitude. Therefore, we expect fn+1 to be quite large, something like 1452/.0002 = 7 million . I have not gotten these numbers calculated correctly, but the point is roughly that fn+1 / fn ~ 1/[!Syntax Error, Idx φn+1(x) Pn+1(x)] . Because this overlap is so small, the fn+1 numbers keep increasing dramatically relative to the previous fn. Final Conclusion: Solve these "integral equations" for f(x): !Syntax Error, Idx f(x) P0(x) = 1 !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4... The solution f(x) is a horrible looking distribution which I think is equivalent to 0, similar in nature to F(x) = limN→∞ Ncos(Nx). So in essence, f(x) = 0 is the only solution to these equations, despite the fact that this seems to conflict with the requirement that area = 1 from the first equation. Certainly if we looked for a solution of the form f(x) = constant, we find there is no solution since the first item requires constant ≠ 0 while the second requires constant = 0, since, for example, !Syntax Error, Idx P3(x) = -1/8 7. My Problem, Plan D (numeric subdivision of the interval) Solve these "integral equations" for f(x): !Syntax Error, Idx f(x) P0(x) = 1 !Syntax Error, Idx f(x) Pn(x) = 0 n = 1,2,3,4... Break up the interval 0 to 1 into N equal pieces of width Δx = 1/N. Then we have Σi=1N f(xi)Pn(xi) Δx = δn,0 For Maple's sake, rewrite this as Σi=1N f(xi)Pn-1(xi) /N = δn,1 Now define xi = (i/N) so that x1 = 1/N and xN = 1 fi = f(xi) Fni = Pn-1(xi)/N Then our equation is this Σi=1N Fni fi = δn,1 which in matrix form is this F f = e1 and f = F-1(e1) which is similar to what we found above in a different context. Here is my code and then some results: > restart; > with(orthopoly); > with(linalg); > N := 10: > FF := (n,i) -> P(n-1,i/N)/N: > F := matrix(N,N,FF): > FI := inverse(F): > with(plots):listplot(col(FI,1)); N = 10 N = 20 N = 30 N = 40 N = 60 (scaled 1e-38) I think this method is also telling me that f(x) = 0 in the limit N=∞ . These "sines" have only one point on each side. You see the burst moving toward the center leaving f ≈ 0 outside the burst. You see the absurdly large values of f(xi). You have to keep increasing Digits to get clean results, up to 50 here. I cannot really do any larger cases, but I can guess what happens. The burst zeros in on the center, gets even larger, maintains that shape I show. This burst is NOT a delta function or any derivative thereof. Infinite crossings, infinite amplitude, infinitely small range at the half way point. In this little model, as I increase the number of points between 0 and 1, I also increase the number of "conditions" on f(x). We need 10 conditions on 10 points in order to have a well-defined solution. Notice that the basis functions method gives a somewhat cleaner view of this problem for a finite N.