# Solve a symbolic underdetermined Linear System

Dear StackExchange Community,

I'm trying to solve an indeterminate linear system of equations, with $n+1$ variables and $n$ equations; therefore, I need to express all $n$ other variables a function of one of them, in the case bellow p[0,0]. Basically I evaluate this system for some numeric values of λ2, μ2, c2, α[·,·], β[·] and a >= 1:

Table[Power[Sqrt[λ2/(μ2 c2)], n] p[0, n + 1] ==
β[n] p[0, 0] +   Sum[α[n, k] p[0, k], {k, 0, a}], {n, 0, a - 1}]


which I'd like to solve for:

Table[p[0, n], {n, 1, a}]


for the parameter p[0,0], in the form p[0,i] = Ci p[0,0], where Ci's are constants.

Currently, I'm doing:

NumVal = {λ1 -> 1.2, μ1 -> 1.0, c1 -> 1.0, λ2 -> 12.0, μ2 -> 1.0, c2 -> 10.0};

plist = Flatten[Solve[Table[
Power[Sqrt[λ2/(μ2 c2)], n] p[0, n+1] ==
β[n] p[0, 0] + Sum[α[n, k] p[0, k], {k, 0, a}] /.
NumVal, {n, 0, a-1}], Table[p[0, n], {n, 1, a}]], 1];


The values of α[·,·], β[·] come from somewhere else.

It was working for small values of a, but now that I want to evaluate it for a=100 Mathematica returns an error. So, I'd like to know if you guys would know a better away for me to compute it which could surpass this limitation.

Thank you very much.

Best,

GT

EDIT 1

Dear,

Still, I was trying to see if I could understand what was going on, but even for small example below Mathematica returns: Solve::svars: Equations may not give solutions for all "solve" variables!. It works up to a = 32 only.

NumVal = {λ1 -> 1.2, μ1 -> 1.0, c1 -> 1.0, λ2 -> 12.0, μ2 -> 1.0, c2 -> 10.0}
rad = Sqrt[λ2/(μ2 c2)] /. NumVal;

atable = {50, 100};
i3 = 1;

For[in = 0, in <= atable[[i3]] - 1, in++, {
β[in] = RandomReal[];
For[ik = 0, ik <= atable[[i3]], ik++, {
α[in, ik] = RandomReal[];
}]
}]

plist = Flatten[
Solve[
Table[
Power[rad, n] p[0, n + 1] ==
β[n] p[0, 0] + Sum[α[n, k] p[0, k], {k, 0, atable[[i3]]}],
{n, 0, atable[[i3]] - 1}],
Table[p[0, n], {n, 1, atable[[i3]]}]
], 1]


What would be the best way to solve this linear system?

Thank you,

GT

• Your plist definition does not seem syntactically correct. Maybe something was missed in copy/pasting. Could you double-check it? Commented Aug 11, 2015 at 15:35
• Dear, @MarcoB, I checked it. Thank you. Commented Aug 11, 2015 at 15:42
• "Mathematica returns an error." What is the error? Commented Aug 11, 2015 at 16:07
• When I export it, I get this: λ1 μ1 c1 λ2 μ2 c2 a ℒ1 ℒ2 Case 1.2 1. 1 12. 1. 10. 0 0.035653 0.179767 2 1.2 1. 1 12. 1. 10. 10 0.075221 0.175811 1 1.2 1. 1 12. 1. 10. 20 0. - 2.837882191385269*^6*p[0, 1] + 4.281087859176636*^6*p[0, 2] - 2.8133612817382812*^6*p[0, 3] + 676003.4604458809*p[0, 4] + (638312.707710743*(2.1999999999999993 + 3.405458629662323*^6*p[0, 1] - (...) 1 Commented Aug 12, 2015 at 7:43
• Isn't that the error (it's actually a warning) that you should expect in an underdetermined system? It will solve for $n$ variables in terms of the remaining variable. I wonder why you don't always get it.... Commented Aug 12, 2015 at 16:49

We will go by solving one equation at a time and generating the corresponding replacement rules. Beware of possible numerical instabilities.

The following is the equivalence between the code in your edited example and my code on the previous incarnation of this answer. I think this is enough for you to use it. Please note that the only claim on the following code is that it seems to work. I believe there should be better ways to do it.

SeedRandom[42];
NumVal = {λ1 -> 1.2, μ1 -> 1.0,  c1 -> 1.0, λ2 -> 12.0, μ2 -> 1.0, c2 -> 10.0};
rad = Sqrt[λ2/(μ2 c2)] /. NumVal;
atable = {10, 31, 32, 50, 100};
i3 = 1;
For[in = 0, in <= atable[[i3]] - 1, in++, {β[in] = RandomReal[];
For[ik = 0, ik <= atable[[i3]],  ik++, {α[in, ik] = RandomReal[];}]}]

yourT = Table[Power[rad, n] p[0, n + 1] == β[n] p[0, 0] +
Sum[α[n, k] p[0, k], {k, 0, atable[[i3]]}], {n, 0, atable[[i3]] - 1}];

plist = Flatten[Solve[yourT, Table[p[0, n], {n, 1, atable[[i3]]}]], 1]
Table[p[0, n], {n, 1, atable[[i3]]}]], 1]
(*
{p[0, 1] -> -0.508959 p[0, 0],  p[0, 2] -> -0.423634 p[0, 0],
p[0, 3] -> -0.745538 p[0, 0],  p[0, 4] -> -1.0122 p[0, 0],
p[0, 5] -> 0.00348381 p[0, 0], p[0, 6] -> -0.057136 p[0, 0],
p[0, 7] -> -0.460527 p[0, 0],  p[0, 8] -> 0.252575 p[0, 0],
p[0, 9] -> 0.259387 p[0, 0],   p[0, 10] -> -0.730669 p[0, 0]}
*)


"My" code:

rul = {p[0, 0] -> p[0, 0]};
Monitor[
For[i = 1, i <= atable[[i3]], i++,
AppendTo[rul, Solve[yourT[[i]] //. rul, p[0, i]][[1, 1]]];
rul = Thread[rul[[All, 1]] -> (rul[[All, 2]] //. rul // Simplify)]], i];

rul

(*
{p[0, 0] -> p[0, 0],
p[0, 1] -> -0.508959 p[0, 0],  p[0, 2]  -> -0.423634 p[0, 0],
p[0, 3] -> -0.745538 p[0, 0],  p[0, 4]  -> -1.0122 p[0, 0],
p[0, 5] -> 0.00348381 p[0, 0], p[0, 6]  -> -0.057136 p[0, 0],
p[0, 7] -> -0.460527 p[0, 0],  p[0, 8]  -> 0.252575 p[0, 0],
p[0, 9] -> 0.259387 p[0, 0],   p[0, 10] -> -0.730669 p[0, 0]}
*)


You may see that the results are the same, but my "method" (solving one equation at a time) can (slowly) process more than 100 variables without much problem.

• Dear, @belisarius, I have a code that is "fragile", and my coding skills in Mathematica are not that high. I could not really understand your answer to use it in my file. I send you my [file] (goo.gl/D8T4LG ), so if eventually you have the time and the endeavour to help me implement your solution. Anyhow, thank you very much for you assistance. Commented Aug 12, 2015 at 9:02
• @GuilhermeThompson See edit please Commented Aug 12, 2015 at 16:22
• Are you sure those For loops won't break my Mathematica? Commented Aug 12, 2015 at 16:57
• @MichaelE2 I believe this is the best way for the OP to understand the solution :D . A mapped/functional approach won't differ too much in this case. Commented Aug 12, 2015 at 17:01

Second update -- I should state in simplest terms the issue the OP is facing.

The set-up. The equations for a given a are

 eqs = Table[Power[Sqrt[λ2/(μ2 c2)], n] p[0, n + 1] == β[n] p[0, 0] +
Sum[α[n, k] p[0, k], {k, 0, a}] /. NumVal, {n, 0, a - 1}];


All the variables involved in eqs are given by

 vars = Table[p[0, n], {n, 0, a}];     (* starts at n == 0 *)


In the OP's code, the variables passed to Solve are

 Table[p[0, n], {n, 1, a}]             (* starts at n == 1, a == atable[[i3]] *)


which is the same as Rest[vars].

The issue. The system eqs is linear with entirely numeric (machine-real) coefficients and so should be solvable quickly. Indeed, even for a = 100, Solve[eqs, vars] is quite fast (< 0.06 sec.); it can be made faster still if solved as a matrix problem. However, in the call Solve[eqs, Rest[vars]] that is equivalent to the OP's use of Solve, the system is now linear with a symbolic parameter, namely p[0, 0]. This slows things down and in the case where a = 100, by quite a lot longer than I'm willing to wait.

Upshot: Requesting Solve to solve for all variables is faster than asking it to solve for exactly the "right number" of variables. No doubt this is because the problem can be formulated in terms of numeric matrices.

The non-issue. When Mathematica emit messages, some are warnings and some are errors. The user has to judge whether a warning is an error or not. Solve[eqs, vars] emits the message

Solve::svars: Equations may not give solutions for all "solve" variables. >>

Since there is one more variable than equations, it is to be expected. Indeed, I take it as a reassurance: to not get the message would indicate an error in the code.

Example. With the set up of belisarius/OP (SeedRandom[42]; a = atable[[5]], i.e., a = 100):

Solve[eqs, vars] // AbsoluteTiming


Solve::svars: Equations may not give solutions for all "solve" variables. >>

(*
{0.056994,
{{p[0, 1] -> 0. - 0.683237 p[0, 0], p[0, 2] -> 0. - 0.180283 p[0, 0],
...,
p[0, 99] -> 0. - 0.0000293422 p[0, 0], p[0, 100] -> 0. + 0.0000153311 p[0, 0]}}}
*)


First update and original answer -- see edit history.