10

I am working with functions calculated from a set of general basis functions.

f = a[x] + c1*b[x] + c2*c[x];
g = Expand[f*f];
h = Integrate[g, {x, -Infinity, Infinity}];

This much works fine, and I get a nice expression with products of combinations of a[x], b[x], and c[x] in the integrand. But now I need to calculate the values of c1 and c2 that optimize h and the optimal value of h. This fails:

cs = {c1, c2};
s = Solve[{D[h, #] & /@ cs == 0}, cs]
(* Solve::nsmet: This system cannot be solved with the methods available to Solve. >> *)

Is there any hack to make this work?

EVEN better would be if I could get the answer in tidy bracket notation, where for example denotes the definite integral of a[x]*b[x].

dionys
  • 4,321
  • 1
  • 19
  • 46
Jerry Guern
  • 4,602
  • 18
  • 47

3 Answers3

16

Similar idea to belisarius, except in V10 we can inactivate Integrate to keep it from evaluating or even trying to evaluate:

h = Inactive[Integrate][g, {x, -Infinity, Infinity}]

It is not necessary in this example, as belisarius' answer shows, but one of its intended uses is to do algebra/calculus on integrals and derivatives. Inactive can be removed easily with

Activate[h]

The function linearExpand expands its argument according to linearity properties. Factors/terms that do not depend on x are treated as constants (see update below for a more general approach).

Clear[linearExpand];
linearExpand[e_] := e //. {int : Inactive[Integrate][_Plus, _] :> Distribute[int], 
    Inactive[Integrate][integrand_Times, dom : {x_, _, _} | x_] :> 
     With[{dependencies = Internal`DependsOnQ[#, x] & /@ List @@ integrand},
      Pick[integrand, dependencies, False] *
       Inactive[Integrate][Pick[integrand, dependencies, True], dom]
      ]};

OP's sample problem:

Solve[D[h, #] == 0 & /@ cs // linearExpand, cs]

Mathematica graphics

D[h, #] == 0 & /@ cs // linearExpand

Mathematica graphics


For what it's worth...

...here's a general linearity expander. Considers factors that do not depend on x, which may be a list of symbols, as constants.

linearExpand[e_, x_, head_] := 
  e //. {op : head[arg_Plus, __] :> Distribute[op], 
    head[arg1_Times, rest__] :> 
     With[{dependencies = Internal`DependsOnQ[#, x] & /@ List @@ arg1},
      Pick[arg1, dependencies, False] head[
        Pick[arg1, dependencies, True], rest]
      ]};

Examples:

linearExpand[D[h, #] == 0 & /@ cs, x, Inactive[Integrate]]
(* same as above *)

linearExpand[foo[(a[x] + c b[y]) (2 a[x] - c b[y]) // Expand, randomarg], x, foo] (* -c^2 b[y]^2 foo[1, randomarg] + c b[y] foo[a[x], randomarg] + 2 foo[a[x]^2, randomarg] *)

linearExpand[foo[(a[x] + c b[y]) (2 a[x] - c b[y]) // Expand, randomarg], {x, y}, foo] (* 2 foo[a[x]^2, randomarg] + c foo[a[x] b[y], randomarg] - c^2 foo[b[y]^2, randomarg] *)

Michael E2
  • 235,386
  • 17
  • 334
  • 747
13
f = a[x] + c1*b[x] + c2*c[x];
g = Expand[f*f];

(* we need to get the constants out of the integrals first*)

h = Distribute@Integrate[g, {x, -∞, ∞}] //. Integrate[q1___  r__  q2___, {v_, s__}] /; 
                                        FreeQ[{r}, v] :> r Integrate[q1 q2, {v, s}];

s= Solve[And @@ Thread[D[h, #] & /@ {c1, c2} == 0], {c1, c2}]

(* now we go to bra-ket notation *)

bkRulez = {Integrate[a_ [x] b_[x], {x, -∞, ∞}]     ->   AngleBracket[a, b], 
           Integrate[Power[a_ [x], 2], {x, -∞, ∞}] ->   AngleBracket[a, a]}
Column @@ (s /. bkRulez) // TeXForm

$$\begin{array}{l} \text{c1}\to -\frac{\langle a,c\rangle \langle b,c\rangle -\langle c,c\rangle \langle a,b\rangle }{\langle b,c\rangle ^2-\langle b,b\rangle \langle c,c\rangle } \\ \text{c2}\to -\frac{\langle b,b\rangle \langle a,c\rangle -\langle a,b\rangle \langle b,c\rangle }{\langle b,b\rangle \langle c,c\rangle -\langle b,c\rangle ^2} \\ \end{array}$$

Dr. belisarius
  • 115,881
  • 13
  • 203
  • 453
  • Cool, thanks! Can you tell me how you got the ouput to format that way on mathematica.SE? When I try to cut-paste output, it comes out as hard-to-read flat text, not formatted like in my mathematica window. The editing-help page strangely does not cover this topic. – Jerry Guern Oct 30 '14 at 17:52
  • @JerryGuern Try enclosing the TeXForm in $ ... $ or in $$ ... $$ – Dr. belisarius Oct 30 '14 at 18:03
  • @JerryGuern Check this http://meta.math.stackexchange.com/questions/5020/mathjax-basic-tutorial-and-quick-reference – Dr. belisarius Oct 30 '14 at 18:07
  • All that did was put the text in a different font. But it's still writing out "Integral" etc instead of putting a nice integral symbol like in the text I copied. – Jerry Guern Oct 30 '14 at 18:34
  • @JerryGuern But you need to put the TeXForm[ expr ] inside the $ ... $ ! – Dr. belisarius Oct 30 '14 at 18:43
  • @JerryGuern See the source of my post here http://mathematica.stackexchange.com/revisions/a6368cdf-7c39-4b1b-9273-94ce840f2d3f/view-source – Dr. belisarius Oct 30 '14 at 18:45
  • Thanks for all the help! But I have a minor coding question, because I am new to Mathematica. Why did you do this: Solve[And @@ Thread[D[h, #] & /@ {c1, c2} == 0], {c1, c2}]
    instead of this?: Solve[{D[h, #] & /@ {c1, c2} == 0}, {c1, c2}] Do the And and Thread buy something?
    – Jerry Guern Oct 30 '14 at 22:18
  • @JerryGuern The form I used here is useful in more general situations where you have AND/OR clauses. For this particular problem your way is enough – Dr. belisarius Oct 30 '14 at 22:31
1

I played with the commands suggested above and found that Distributed[] solved the first half of the problem but only for indefinite integrals:

f = a[x] + c1*b[x] + c2*c[x];
g = Expand[f*f];
h = Distribute[Integrate[g, x]];
cs = {c1, c2};
Solve[{D[h, #] & /@ cs == 0}, cs]

I worry when my solution looks too simple... Are the more complicated answers above handling some problem I'm not aware of?

Jerry Guern
  • 4,602
  • 18
  • 47
  • Your solution doesn't seem to work with your original equation Solve[{D[h, #] & /@ {c1, c2} == 0}, {c1, c2}], which is necessary if you plan to have more cis... – Dr. belisarius Oct 30 '14 at 18:27
  • @belisarius Why did you say it doesn't work? It worked fine for me. But I edited the code above so the Solve line can handle a long cs list. – Jerry Guern Oct 30 '14 at 18:45
  • 1
    Well ... it didn't work for me yesterday :( Or so I believe – Dr. belisarius Oct 30 '14 at 18:47
  • @JerryGuern Neat. I'm sure I tried that and the constants did not come out of the integrals (for h). Then Solve wouldn't work. – Michael E2 Oct 30 '14 at 18:47
  • @belisarius It's an odd/even bug. Works on even-numbered days. – Michael E2 Oct 30 '14 at 18:48
  • 1
    @MichaelE2 No, I got it. Now Jerry is working with indefinite integrals and it works. Not for definite integrals, though – Dr. belisarius Oct 30 '14 at 18:50
  • @belisarius Dang it, you're right, it doesn't work for definite integrals. Well, this is why I keep it humble and tentative when I question experts. – Jerry Guern Oct 30 '14 at 19:08
  • @belisarius If the domain of integration is fixed throughout the problem, Jerry could probably do the algebra with indefinite integrals. They seem to work better. – Michael E2 Oct 30 '14 at 19:37
  • @belisarius One might add such a construct h=Distribute[Integrate[g, x]] /. Integrate[k_, x] -> HoldForm[Integrate[k, {x, -Infinity, Infinity}]] // ReleaseHold and it will enable one to work with definite integrals. But the question is WHY it does not work for definite integrals plainly? – Alexei Boulbitch Aug 14 '15 at 08:33