Prev: Using Mathematica for Electrophysiology Data Analysis
Next: formatted table output to ascii file
From: leigh pascoe on 16 Jun 2010 05:41
Dear Experts/Adam Strebonski,

I have been trying for several years to solve a set of simultaneous
polynomial equations. Since the algorithms used by Mathematica are
guaranteed to find solutions to these types of equations it should be a
slam-dunk. However I have been frustrated by memory limitations (Windows
XP/quad core with 4 Go) - Mathematica crashes out with a memory exceeded
error.

Having seen the post "Big Memory needed" and the welcome development
that V8 of Mathematica is much better at these problems, I am hoping
that someone with access to V8 and a large capacity machine may be able
to help me. The equations are given below and occur in the theory of
population genetics. The variables of interest are denoted x1, x2, x3, x4

First define the variables

d := x2*x3 - x1*x4
c := r*w14*d

and some constants (symmetric matrix)

w21 := w12
w31 := w13
w41 := w14
w23 := w14
w42 := w24
w43 := w34
w32 := w14

We use the following intermediate functions to simplify the writing of
the equations:

w1 := x1*w11 + x2*w12 + x3*w13 + x4*w14
w2 := x1*w12 + x2*w22 + x3*w23 + x4*w24
w3 := x1*w13 + x2*w23 + x3*w33 + x4*w34
w4 := x1*w14 + x2*w24 + x3*w34 + x4*w44

Now the equations can be written

eq1 = x1 + x2 + x3 + x4 == 1
eq2 = x1*(x1*w1 + x2*w2 + x3*w3 + x4*w4) == x1*w1 + c
eq3 = x2*(x1*w1 + x2*w2 + x3*w3 + x4*w4) == x2*w2 - c
eq4 = x3*(x1*w1 + x2*w2 + x3*w3 + x4*w4) == x3*w3 - c

These are three equations of 3rd degree in x's and a linear constraint.
I attempt to solve them with the Solve command

FullSimplify[Solve[{eq1, eq2, eq3, eq4}, {x1, x2, x3, x4}]]

No more memory available.
Mathematica kernel has shut down.
Try quitting other applications and then retry.

If anyone can help me to know the number of solutions and their
formulae, I would be greatly appreciative. Eventually I would also like
to subject the solutions to the restriction that all the xi's are in the
closed interval [0,1], which of course depends on the values of the
constants.

Mathematica is able to find the solutions of the special case when c=0

eq5 = x1 + x2 + x3 + x4 == 1
eq6 = x1*(x1*w1 + x2*w2 + x3*w3 + x4*w4) == x1*w1
eq7 = x2*(x1*w1 + x2*w2 + x3*w3 + x4*w4) == x2*w2
eq8 = x3*(x1*w1 + x2*w2 + x3*w3 + x4*w4) == x3*w3

Timing[FullSimplify[Solve[{eq5, eq6, eq7, eq8}, {x1, x2, x3, x4}]]]

The 15 solutions to the simplified case are found in about 8 secs and
occupy several screens. However these solutions can be found by hand and
written much more succintly than the Mathematica result. Thanks in
advance for any help.

LP

 | 
Pages: 1
Prev: Using Mathematica for Electrophysiology Data Analysis
Next: formatted table output to ascii file