How do I solve N simultaneous equations for N variables?

Question

I have a function:

f[x_] := x + 31 x^3 + 5 x^25

Which I want to find an expansion for:

ex[x_, a_, b_, c_, d_] := a + b x + c x^3 + d x^25

I do this by solving to find the coefficients a,b,c,d by setting the expansion equal to the function at the appropriate number of points:

sol = Solve[{ex[1, a, b, c, d] == f[1], ex[3, a, b, c, d] == f[3], 
ex[2, a, b, c, d] == f[2], ex[4, a, b, c, d] == f[4]}]

This works fine, but its ugly. I've only used this to illustrate what I want to do... but actually I want to find an expansion involving 'N' expansion coefficients, and perform solve at N different points to find what they are. However, I don't understand how to do that without having functions with N+1 arguments and so on.

If I haven't expressed myself clearly let me know and I'll add more.

EDIT: Perhaps I'll add that really the expansion I'm looking for is something like (though how I put the coefficients into the function argument I don't know):

ex[x_,a[i]_,b[i]_,N_,M_]:= Sum[a[n] ChebyshevT[n,x],{n,0,N}] + Sum[b[n] * y[n,x],{n,1,M}]

Where the functions y[n,x] are some other functions which is not orthogonal to the Chebyshev polynomials.

Answer 1

(Edited; original at the end)

Gram-Schmidt orthogonalization provides an answer. Let's use a running example to illustrate. It begins even before the Chebyshev polynomials, with their domain--the interval $[-1,1]$--and the kernel for which they are orthogonal:

limits = {-1, 1};
k[x_] := 2/Sqrt[1 - x^2]/Pi;

The Chebyshev polynomials are obtained by means of the Gram-Schmidt orthogonalization process from the monomial functions $(1, x, x^2, \ldots, x^n, \ldots)$. For instance, here are the first four of them as obtained from the powers $0$ through $3$:

 Orthogonalize[{1, x, x^2, x^3}, Integrate[#1 #2 k[x], {x, limits[[1]], limits[[2]]}] &]  
  // Expand

$\left\{\frac{1}{\sqrt{2}},x,2 x^2-1,4 x^3-3 x\right\}$

(If you compare this to ChebyshevT[#, x] & /@ Range[0, 3] // TraditionalForm you will see that the first is rescaled but all the rest are the same; this is because Mathematica for some reason does not use an orthonormal basis for ChebyshevT.)

Thus, in order to obtain an expansion in terms of a finite set of Chebyshev polynomials and other functions, first orthogonalize them, normalize those if necessary, and obtain the coefficients in terms of the inner products.

There's a problem: unless the collection of functions is very nice, Mathematica will be unable to obtain closed-form expressions for the orthogonal basis. Instead, we can numerically integrate.

For example, let's start with a few power functions--they will generate Chebyshev polynomials automatically--together with some stranger transcendental ones (plotted below):

vectors = Function /@ Table[#^k, {k, 0, 6}];
y = {Sqrt[1 - Abs[#]] &, Exp, # Abs[Log[#]] &};
Plot[Evaluate@Through[y[x]], {x, limits[[1]], limits[[2]]}]

In this fashion we have recreated the situation presented in the problem: we seek to expand arbitrary functions in terms of Chebyshev polynomials--which are mutually orthogonal--and some other functions y which are not necessarily orthogonal to the Chebyshev polynomials.

Obtaining an orthogonal basis takes a few seconds and some fiddling to compute the integrals with sufficient accuracy:

t = Orthogonalize[Through[(vectors~Join~y)[x]], 
 NIntegrate[#1 #2 k[x], {x, limits[[1]], limits[[2]]}, 
   MaxRecursion -> 12] &] // Simplify // Chop // Rationalize

Let's check that this basis is really orthonormal (if not, we would have to normalize its elements):

Outer[NIntegrate[#1 #2 k[x], {x, limits[[1]], limits[[2]]}, MaxRecursion -> 12] &, t, t]

I applied Chop[%, 10^-6] to remove some small near-zero values and obtained

$$\left( \begin{array}{cccccccccc} 1. & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1. & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1. & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1. & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1. & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1. & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1. & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1. & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1. \end{array} \right)$$

The off-diagonal zeros confirm orthogonality. The lengths of the vectors are on the diagonal: already normalized, as intended. Notice the zero on the diagonal! To within the limits of numerical error, it appears one of these functions is a linear combination of the others. Nevertheless, we can still proceed without having to take any special steps. Let's pause first to look at these basis functions:

GraphicsGrid[{Plot[Evaluate@#, {x, limits[[1]], limits[[2]]}] & /@ Partition[t, 5]}]

Definitely some of the later ones (on the right) aren't pure Chebyshev polynomials!

Let's expand a function in terms of this basis; some function that clearly is not within the linear span of the original set of functions:

f[x_] := x + 31 x^3 + 5 x^25
a = NIntegrate[f[x] # k[x], {x, limits[[1]], limits[[2]]}] & /@ t // Chop;
a . t // Simplify

$2.29258 \sqrt{1-|x|}-6.43222 x |\log (x)|-3.82443 x^6+21.3062 x^5+10.7913 x^4+6.02688 x^3-14.6504 x^2+18.7429 x-2.6231$

(Numerical error is not as critical in this calculation, so we don't have to be so fiddly.) The vector a is the desired set of expansion coefficients. Its inner product with t expresses the expanded value: it is the orthogonal projection of $f$ onto the space of functions spanned by the functions we defined earlier in vectors and y. ("Orthogonal" means relative to the kernel $k$.)

How does the original function compare to its expansion?

Plot[{f[x], a.t}, {x, limits[[1]], limits[[2]]}]

This looks pretty good. However, because f gets large, we can't see small differences. Let's look at the residuals of the expansion instead. There are two meaningful ways to do this: (a) by subtracting f from its expansion and (b) by dividing those differences by the values of the kernel: these "standardized" differences ought to be fairly uniform across the domain. Here are both plots:

GraphicsGrid[{{Plot[a.Through[nCn[x]] - f[x], {x, limits[[1]], limits[[2]]}, PlotStyle -> Thick], 
   Plot[ (a.Through[nCn[x]] - f[x]) / k[x], {x, limits[[1]], limits[[2]]}]}}]

These errors are pretty small compared to the original range of $f$ and, as expected, the standardized errors are reasonably uniform.

---

Original Reply

Chebyshev polynomials $\psi_n(x)$ are orthogonal on $[-1,1]$ with respect to $\frac{dx}{\sqrt{1-x^2}}$. The normalization used by Mathematica makes all but the zeroth have norm $\pi/2$; the zeroth has norm $\pi$; e.g.,

Integrate[ChebyshevT[#, x]^2 / Sqrt[1 - x^2] , {x, -1, 1}] & /@ Range[0, 5]

$\left\{\pi ,\frac{\pi }{2},\frac{\pi }{2},\frac{\pi }{2},\frac{\pi }{2},\frac{\pi }{2}\right\}$

Brute Force Method

Therefore the coefficient of the $n$th Chebyshev polynomial is found by integrating $f$ against ChebyshevT[n, x] / Sqrt[1 - x^2] and dividing by $\pi/2$ (by $\pi$ when $n=0$). (For an odd function $f$ only the odd coefficients will be nonzero: knowing this can halve the execution time.) Only the Chebyshev polynomials up to the degree of $f$ have to be evaluated:

f[x_] := x + 31 x^3 + 5 x^25;
degree = Length[CoefficientList[p[x], x]] - 1;
a = Integrate[(2/Pi) p[x] ChebyshevT[#, x] / Sqrt[1 - x^2] , {x, -1, 1}] & /@ Range[0, degree];
a = a Prepend[ConstantArray[1, degree], 1/2]

$\left\{0,\frac{108212247}{4194304},0,\frac{19038803}{2097152},0,\frac{2042975}{2097152}, \ldots, \frac{125}{16777216},0,\frac{5}{16777216}\right\}$

As a check, compute the expansion to verify it equals the original polynomial:

a . (ChebyshevT[#, x] & /@ Range[0, degree]) // Expand

$5 x^{25}+31 x^3+x$

Fast Method

Using the preceding, it's straightforward to check and prove that the coefficient of the $n$th Chebyshev polynomial in the expansion of $x^k$ is related to a Binomial coefficient:

Clear[c];
c[k_, n_] /; n <= k && EvenQ[n - k] := c[k,n] 
  = 2^(1 - k) Binomial[k, (k - n)/2] If[n == 0, 1/2, 1];
c[k_, n_] := 0

Whence, to expand any polynomial f, just replace the powers of its argument by their expansions:

f[x] /. {Times[a_, Power[x, e_]] :> 
 a Sum[c[e, k] Subscript[\[Psi], k][x], {k, Mod[e, 2], e, 2}], 
 x -> Subscript[\[Psi], 1][x]} // Expand 

$\frac{108212247 \psi _1(x)}{4194304}+\frac{19038803 \psi _3(x)}{2097152}+\frac{2042975 \psi _5(x)}{2097152}+\ldots+\frac{125 \psi _{23}(x)}{16777216}+\frac{5 \psi _{25}(x)}{16777216}$

---

A similar approach will work with other systems of orthogonal polynomials.

Answer 2

You could do something like:

b[i_, x_] := ChebyshevT[i - 1, x];(*Tch basis as an example*)
app[f_, pSet_]:= 
     Solve[Table[Sum[a@i b[i, x], {i,Length@#}]==f@x, {x, #}], Array[a,Length@#]] &@  pSet

Where pSet is the list of points (x-values) you use for equaling the functions.

Use it like this:

app[# + 31 #^3 + 5 #^25 &, Range@5]

Edit

Note that for example:

app[# + 31 #^3 + 5 #^25 &, Range@3]

Expands to:

Solve[{ a[1] +   a[2] +    a[3] == 37, 
        a[1] + 2 a[2] + 7  a[3] == 167772410, 
        a[1] + 3 a[2] + 17 a[3] == 4236443048055}, 
     {a[1], a[2], a[3]}]

So you could also try some linear method directly