# How do I solve N simultaneous equations for N variables?

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.

-
Is your function always a polynomial ? –  b.gatessucks Feb 12 '13 at 15:38
The function isn't but the expansion basis functions will be Chebyshev Polynomials, I just used that to illustrate my point. –  user5866 Feb 12 '13 at 15:40
If your basis isn't orthogonal, your expansion won't (probably) be unique ... –  belisarius Feb 12 '13 at 16:11
As in more than one set of a[n] and b[n] will give f[x] to arbitrary accuracy? (that wouldn't necessarily be a problem but I hadn't realised that) –  user5866 Feb 12 '13 at 16:16
Something like Solve[Table[ex[j, a, b, c, d] == f[j],{j,nn+mm}],Join[Table[a[n],{n,nn}],Table[b[n],{n,mm}]]]? –  Daniel Lichtblau Feb 12 '13 at 16:21

### (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.

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.

-
I posted a similar answer earlier and deleted it because the OP commented Sorry, that is completely true, but it won't work for the actual expansion I want to do which is a mixture of Chebyshev and other singular functions. I'm sorry, i'll add more information, I was just trying to keep it simple for readability –  belisarius Feb 12 '13 at 22:06
@belisarius Thank you for the clarification. Because the techniques in this answer might still have a bearing on the reformulated question, I will leave it up in the meantime. –  whuber Feb 12 '13 at 22:57
But .. why don't you post a question about that and then answer it yourself? The OP's problem is different! –  belisarius Feb 12 '13 at 23:29
@belisarius That's a good suggestion, but I still think that what I originally wrote contributes to the question. Instead, I have--I believe--now fully answered the question in an edit. I have retained the original because for some readers it might help clarify or streamline the application of the first part. –  whuber Feb 13 '13 at 19:25

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

-
Thanks a lot for your response What is the role of pSet? Is it just a list of points where we set the expansion equal to the function? If so can I just use pSet=Table[n,{n,1,N}] to give a list of N points? also should the 1st and 2nd arguments of your Solve function be enclosed in { } or not? –  user5866 Feb 12 '13 at 16:56
never mind I think its a pointer to a set, and you've used Range in the place of my suggested table. so much still to learn –  user5866 Feb 12 '13 at 17:09
@user5866 pSet is a ... points set :) –  belisarius Feb 12 '13 at 17:35
@user5866 There are no pointers natively in Mathematica –  belisarius Feb 12 '13 at 18:25
This works as an answer, but I lack the intelligence to expand it to the case I actually need, which isn't just a 1,2,3 range but rather N set of points. What is it that you are passing with the Range@3 argument? I had it in my head that this was the address of a range of 3 numbers that you were pointing to with pSet.... :S –  user5866 Feb 12 '13 at 19:49