Thread Subject: Solve Equation Symbolically

Good afternoon all,

I am trying to solve the equation cv3 below using the solve
command but for some reason Matlab is not able to solve for
x. Or it returns strange constants in the result. It also
sometimes is dropping some of the symbols. I am just trying
to find out what x equals to. Any help would be greatly
appreciated.
cv0, p[n], x, db and cr are all symbols

cv1 = (cv0 + p1 - (x)*(db-cv0))*(1+cr)
cv2 = (cv1 + p2 - (x)*(db-cv1))*(1+cr)
cv3 = (cv2 + p3 - (x)*(db-cv2))*(1+cr)

I was using solve(cv3, 'x') but this is giving strange
constants.

In article <fha99r$osc$1@fred.mathworks.com>,
evan klinger <eklinger1-nospam@cox.net> wrote:
>Good afternoon all,
>
>I am trying to solve the equation cv3 below using the solve
>command but for some reason Matlab is not able to solve for
>x. Or it returns strange constants in the result. It also
>sometimes is dropping some of the symbols. I am just trying
>to find out what x equals to. Any help would be greatly
>appreciated.
>cv0, p[n], x, db and cr are all symbols

>cv1 = (cv0 + p1 - (x)*(db-cv0))*(1+cr)
>cv2 = (cv1 + p2 - (x)*(db-cv1))*(1+cr)
>cv3 = (cv2 + p3 - (x)*(db-cv2))*(1+cr)

>I was using solve(cv3, 'x') but this is giving strange
>constants.

It's a cubic in x, so there are three solutions. Two of them
involve imaginary numbers. Because of the numbers of terms involved,
each of the three solutions runs to pages and pages.

If you were using maple directly, you would,

> collect(expand(cv3),x);

(-db+3*cv0*cr^2-3*db*cr+3*cv0*cr+cv0*cr^3-3*db*cr^2-db*cr^3+cv0)*x^3+(-2*db*cr
^3+3*cv0-8*db*cr+3*p1*cr^2+p1-7*db*cr^2-3*db+p1*cr^3+3*p1*cr+3*cv0*cr^3+9*cv0*
cr^2+9*cv0*cr)*x^2+(-3*db+3*cv0-4*db*cr^2+9*cv0*cr^2-db*cr^3-6*db*cr+9*cv0*cr+
3*cv0*cr^3+2*p1+6*p1*cr+p2+6*p1*cr^2+2*p2*cr+2*p1*cr^3+p2*cr^2)*x+cv0+p2+p3+p1
+p3*cr+cv0*cr^3+3*cv0*cr+3*p1*cr+3*p1*cr^2+2*p2*cr+p1*cr^3+3*cv0*cr^2+p2*cr^2

Now treat this as a general cubic, A*x^3+B*x^2+C*x+D:

> [solve(A*x^3+B*x^2+C*x+D,x)];
[1/6/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2
*A^2+4*D*B^3)^(1/2)*A)^(1/3)-2/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A^2-8*B^3+12*3^
(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)-1/3*B/A,

-1/12/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^
2*A^2+4*D*B^3)^(1/2)*A)^(1/3)+1/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A^2-8*B^3+12*3
^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)-1/3*B/A+
1/2*I*3^(1/2)*(1/6/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*
C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)+2/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A
^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(
1/3)),

 -1/12/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*
D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)+1/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A^2-8*B
^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)-1
/3*B/A-1/2*I*3^(1/2)*(1/6/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*
B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)+2/3*(3*C*A-B^2)/A/(36*C*B*A-\
108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/
2)*A)^(1/3))]

The second and third of these involve imaginary numbers, so select
the first of them,

1/6/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2
*A^2+4*D*B^3)^(1/2)*A)^(1/3)-2/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A^2-8*B^3+12*3^
(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)-1/3*B/A

and substitute by comparison with the expanded cv3,

A = (-db+3*cv0*cr^2-3*db*cr+3*cv0*cr+cv0*cr^3-3*db*cr^2-db*cr^3+cv0)
B = (-2*db*cr^3+3*cv0-8*db*cr+3*p1*cr^2+p1-7*db*cr^2-3*db+p1*cr^3+3*p1*cr +
     3*cv0*cr^3+9*cv0*cr^2+9*cv0*cr)
C = (-3*db+3*cv0-4*db*cr^2+9*cv0*cr^2-db*cr^3-6*db*cr+9*cv0*cr +
     3*cv0*cr^3+2*p1+6*p1*cr+p2+6*p1*cr^2+2*p2*cr+2*p1*cr^3+p2*cr^2)
D = cv0+p2+p3+p1+p3*cr+cv0*cr^3+3*cv0*cr+3*p1*cr+3*p1*cr^2+2*p2*cr +
    p1*cr^3+3*cv0*cr^2+p2*cr^2

and you will have your multi-page solution.

(Somehow, I suspect that after seeing the multipage solution, you
will decide that maybe solving for x didn't turn out to be
very productive.)
--
   "Okay, buzzwords only. Two syllables, tops." -- Laurie Anderson

I guess I'm confused about why there are 3 solutions since
it can really be written as one long equation if you
substitute the cv's in. In other words, you can write cv3 in
terms of cv0s. Am I missing something?

Thanks!


roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
message <fhadf3$o58$1@canopus.cc.umanitoba.ca>...
> In article <fha99r$osc$1@fred.mathworks.com>,
> evan klinger <eklinger1-nospam@cox.net> wrote:
> >Good afternoon all,
> >
> >I am trying to solve the equation cv3 below using the solve
> >command but for some reason Matlab is not able to solve for
> >x. Or it returns strange constants in the result. It also
> >sometimes is dropping some of the symbols. I am just trying
> >to find out what x equals to. Any help would be greatly
> >appreciated.
> >cv0, p[n], x, db and cr are all symbols
>
> >cv1 = (cv0 + p1 - (x)*(db-cv0))*(1+cr)
> >cv2 = (cv1 + p2 - (x)*(db-cv1))*(1+cr)
> >cv3 = (cv2 + p3 - (x)*(db-cv2))*(1+cr)
>
> >I was using solve(cv3, 'x') but this is giving strange
> >constants.
>
> It's a cubic in x, so there are three solutions. Two of them
> involve imaginary numbers. Because of the numbers of terms
involved,
> each of the three solutions runs to pages and pages.
>
> If you were using maple directly, you would,
>
> > collect(expand(cv3),x);
>
>
(-db+3*cv0*cr^2-3*db*cr+3*cv0*cr+cv0*cr^3-3*db*cr^2-db*cr^3+cv0)*x^3+(-2*db*cr
>
^3+3*cv0-8*db*cr+3*p1*cr^2+p1-7*db*cr^2-3*db+p1*cr^3+3*p1*cr+3*cv0*cr^3+9*cv0*
>
cr^2+9*cv0*cr)*x^2+(-3*db+3*cv0-4*db*cr^2+9*cv0*cr^2-db*cr^3-6*db*cr+9*cv0*cr+
>
3*cv0*cr^3+2*p1+6*p1*cr+p2+6*p1*cr^2+2*p2*cr+2*p1*cr^3+p2*cr^2)*x+cv0+p2+p3+p1
>
+p3*cr+cv0*cr^3+3*cv0*cr+3*p1*cr+3*p1*cr^2+2*p2*cr+p1*cr^3+3*cv0*cr^2+p2*cr^2
>
> Now treat this as a general cubic, A*x^3+B*x^2+C*x+D:
>
> > [solve(A*x^3+B*x^2+C*x+D,x)];
>
[1/6/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2
>
*A^2+4*D*B^3)^(1/2)*A)^(1/3)-2/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A^2-8*B^3+12*3^
>
(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)-1/3*B/A,
>
>
-1/12/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^
>
2*A^2+4*D*B^3)^(1/2)*A)^(1/3)+1/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A^2-8*B^3+12*3
>
^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)-1/3*B/A+
>
1/2*I*3^(1/2)*(1/6/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*
>
C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)+2/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A
>
^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(
> 1/3)),
>
>
-1/12/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*
>
D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)+1/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A^2-8*B
>
^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)-1
>
/3*B/A-1/2*I*3^(1/2)*(1/6/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*
>
B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)+2/3*(3*C*A-B^2)/A/(36*C*B*A-\
>
108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/
> 2)*A)^(1/3))]
>
> The second and third of these involve imaginary numbers,
so select
> the first of them,
>
>
1/6/A*(36*C*B*A-108*D*A^2-8*B^3+12*3^(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2
>
*A^2+4*D*B^3)^(1/2)*A)^(1/3)-2/3*(3*C*A-B^2)/A/(36*C*B*A-108*D*A^2-8*B^3+12*3^
>
(1/2)*(4*C^3*A-C^2*B^2-18*C*B*A*D+27*D^2*A^2+4*D*B^3)^(1/2)*A)^(1/3)-1/3*B/A
>
> and substitute by comparison with the expanded cv3,
>
> A =
(-db+3*cv0*cr^2-3*db*cr+3*cv0*cr+cv0*cr^3-3*db*cr^2-db*cr^3+cv0)
> B =
(-2*db*cr^3+3*cv0-8*db*cr+3*p1*cr^2+p1-7*db*cr^2-3*db+p1*cr^3+3*p1*cr
+
> 3*cv0*cr^3+9*cv0*cr^2+9*cv0*cr)
> C =
(-3*db+3*cv0-4*db*cr^2+9*cv0*cr^2-db*cr^3-6*db*cr+9*cv0*cr +
>
3*cv0*cr^3+2*p1+6*p1*cr+p2+6*p1*cr^2+2*p2*cr+2*p1*cr^3+p2*cr^2)
> D =
cv0+p2+p3+p1+p3*cr+cv0*cr^3+3*cv0*cr+3*p1*cr+3*p1*cr^2+2*p2*cr +
> p1*cr^3+3*cv0*cr^2+p2*cr^2
>
> and you will have your multi-page solution.
>
> (Somehow, I suspect that after seeing the multipage
solution, you
> will decide that maybe solving for x didn't turn out to be
> very productive.)
> --
> "Okay, buzzwords only. Two syllables, tops." -- Laurie
Anderson

"evan klinger" <eklinger1-nospam@cox.net> wrote in message
news:fha99r$osc$1@fred.mathworks.com...
> Good afternoon all,
>
> I am trying to solve the equation cv3 below using the solve
> command but for some reason Matlab is not able to solve for
> x. Or it returns strange constants in the result. It also
> sometimes is dropping some of the symbols. I am just trying
> to find out what x equals to. Any help would be greatly
> appreciated.
> cv0, p[n], x, db and cr are all symbols
>
> cv1 = (cv0 + p1 - (x)*(db-cv0))*(1+cr)
> cv2 = (cv1 + p2 - (x)*(db-cv1))*(1+cr)
> cv3 = (cv2 + p3 - (x)*(db-cv2))*(1+cr)
>
> I was using solve(cv3, 'x') but this is giving strange
> constants.

Since MATLAB uses 'x' as the default independent variable all you need to do
is:

syms cv0 p1 p2 p3 x db cr
cv1 = (cv0 + p1 - (x)*(db-cv0))*(1+cr)
cv2 = (cv1 + p2 - (x)*(db-cv1))*(1+cr)
cv3 = (cv2 + p3 - (x)*(db-cv2))*(1+cr)
solve(collect(expand(cv3))

Hope that helps,
-Ben Hinkle

> solve(collect(expand(cv3))

Sorry for the double post but I tried just solve(cv3) and go the same
solutions so I'm not sure why you are seeing incomplete solutions.

Is there a way to increase the 25,000 character limit? I
really need to solve up to cv12, I just used up to cv3 to
demonstrate what I'm trying to do. However, using
solve(collect(expand(cv12))) returns an error:

??? Error using ==> solve at 162
Unable to find closed form solution.

Error in ==> sym.solve at 49
[varargout{1:max(1,nargout)}] = solve(S{:});

Thanks
Evan

"Ben Hinkle" <bhinkle@mathworks.com> wrote in message
<fhafc8$h1$1@fred.mathworks.com>...
>
> > solve(collect(expand(cv3))
>
> Sorry for the double post but I tried just solve(cv3) and
go the same
> solutions so I'm not sure why you are seeing incomplete
solutions.
>
>

In article <fhaes6$k0h$1@fred.mathworks.com>,
evan klinger <eklinger1-nospam@cox.net> top-posted:

Please do not post your reply at the top of the material you
are commenting on: it makes it more difficult to hold discussions.
Your posting may be the only posting of interest to -you- today,
but I'm following literally over 100 simultaneous discussions today.

>roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
>message <fhadf3$o58$1@canopus.cc.umanitoba.ca>...
>> In article <fha99r$osc$1@fred.mathworks.com>,
>> evan klinger <eklinger1-nospam@cox.net> wrote:

>> >cv0, p[n], x, db and cr are all symbols

>> >cv1 = (cv0 + p1 - (x)*(db-cv0))*(1+cr)
>> >cv2 = (cv1 + p2 - (x)*(db-cv1))*(1+cr)
>> >cv3 = (cv2 + p3 - (x)*(db-cv2))*(1+cr)

>> It's a cubic in x, so there are three solutions. Two of them
>> involve imaginary numbers.

>I guess I'm confused about why there are 3 solutions since
>it can really be written as one long equation if you
>substitute the cv's in. In other words, you can write cv3 in
>terms of cv0s. Am I missing something?

Look at cv1: it has x times cv0.

Now look at cv2: it has x times cv1, so that means it is going
to have x times (x times cv0), or x^2 times cv0.

Now look at cv3: it has x times cv2, so that means it is going
to have x times (x^2 times cv0), or x^3 times cv0.

Therefor your cv3 is going to involve x^3, and is thus
a cubic polynomial. When you solve a cubic polynomial for 0,
then according to the Fundamental Theorem of Algebra, you
are going to have 3 roots of the polynomial.

If your constants -happen- to have very specific (and complicated)
relationships, then all three of the roots would be real, so
for completeness you really shouldn't ignore the two solutions
that involve imaginary numbers: those are of the form
f(x) + g(x)*I and if g(x) just happens to come out to 0 then
you have another real solution.
--
  "I will speculate that [...] applications [...] could actually see a
  performance boost for most users by going dual-core [...] because it
  is running the adware and spyware that [...] are otherwise slowing
  down the single CPU that user has today" -- Herb Sutter

In article <fhag57$bca$1@fred.mathworks.com>,
evan klinger <eklinger1-nospam@cox.net> wrote:
>Is there a way to increase the 25,000 character limit? I
>really need to solve up to cv12, I just used up to cv3 to
>demonstrate what I'm trying to do. However, using
>solve(collect(expand(cv12))) returns an error:

>??? Error using ==> solve at 162
>Unable to find closed form solution.

It isn't any 25000 character limit that is getting you:
you are trying to find the closed form solution to a polynomial
of order 12 that does not happen to be factorable down
to polynomials of degree 5 or less (degree 5 is the maximum
for which closed form polynomial solutions always exist.)

Using Maple directly, there is simply no closed form solution
available for cv12. Beyond cv3, even where the closed form
solutions exist in principle, Maple will just show a generalized
placeholder instead (closed form quadratics and quintics are
pretty messy at the best of times.)
--
  "Is there any thing whereof it may be said, See, this is new? It hath
  been already of old time, which was before us." -- Ecclesiastes

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
message <fhaiam$k8$1@canopus.cc.umanitoba.ca>...
> In article <fhag57$bca$1@fred.mathworks.com>,
> evan klinger <eklinger1-nospam@cox.net> wrote:
> >Is there a way to increase the 25,000 character limit? I
> >really need to solve up to cv12, I just used up to cv3 to
> >demonstrate what I'm trying to do. However, using
> >solve(collect(expand(cv12))) returns an error:
>
> >??? Error using ==> solve at 162
> >Unable to find closed form solution.
>
> It isn't any 25000 character limit that is getting you:
> you are trying to find the closed form solution to a
polynomial
> of order 12 that does not happen to be factorable down
> to polynomials of degree 5 or less (degree 5 is the maximum
> for which closed form polynomial solutions always exist.)
>
> Using Maple directly, there is simply no closed form solution
> available for cv12. Beyond cv3, even where the closed form
> solutions exist in principle, Maple will just show a
generalized
> placeholder instead (closed form quadratics and quintics are
> pretty messy at the best of times.)
> --
> "Is there any thing whereof it may be said, See, this is
new? It hath
> been already of old time, which was before us." --
Ecclesiastes

So what are my options here? I would like a closed form
solution so that I can plug in values into the variables
directly and solve, but if that's not possible, is there any
alternative?

Thank you for all your help so far

In article <fhajsd$bcg$1@fred.mathworks.com>,
evan klinger <eklinger1-nospam@cox.net> wrote:
>roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
>message <fhaiam$k8$1@canopus.cc.umanitoba.ca>...
>> In article <fhag57$bca$1@fred.mathworks.com>,
>> evan klinger <eklinger1-nospam@cox.net> wrote:
>> >Is there a way to increase the 25,000 character limit? I
>> >really need to solve up to cv12,

>> you are trying to find the closed form solution to a
>polynomial
>> of order 12 that does not happen to be factorable down

>So what are my options here? I would like a closed form
>solution so that I can plug in values into the variables
>directly and solve, but if that's not possible, is there any
>alternative?

What kind of ranges of values can your constants take?
I've been playing around with some test values, and it looks to
me as if you could easily encounter severe problems with loss of
precision in your calculations, especially if cr > 2.

With random values between -1 and 1 for the cv12 variables, it
was fairly common in my trials for there to be an x root somewhere
just a bit more than 2, but that by x = 3 that the function value
had risen beyond 10^8 or 10^9 -- but with another one of my trials,
the roots ranged from -7 to +4 .

Unless your values (especially cr) are fairly small, I suspect
you will lose precision so quickly as to make it very doubtful that
you would be able to accurately find the equation roots.

cr is the biggest contribution to the problem: as well as
cv12 being a 12th order polynomial in x, it is also a 12th order
polynomial in cr (cr and x will have the same polynomial order
as each other in each cv*). Thus even a small change in cr can
make a big change in the x roots. And a small change in cr can
result in -different- 12th order roots becoming the real-valued
roots, so even if you had an closed form solution for one of the
12th-order roots, that root might become imaginary with a small
change in parameters.

For example, in a particular trial I did
with random values for the parameters except p3, for p3=1/3,
the first and 7th roots were the real roots, but for p3=1/2,
the 6th and 7th roots were the real roots -- and p3 appears only
in its first power in cv12.


If I recall correctly, I have seen mention of a FEX contribution
that generalizes fzero() to find multiple roots. You are
probably going to have to use something like that, and to
figure out what you want to do with the multiple roots that are
found, but with some kinds of checks in place to ensure that you
are not finding spurious roots (or failing to find valid roots)
due to loss of precision.
--
We regret to announce that sub-millibarn resolution bio-hyperdimensional
plasmatic space polyimaging has been delayed until the release
of Windows Vista SP2.

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
message <fhb2v2$jqi$1@canopus.cc.umanitoba.ca>...
> In article <fhajsd$bcg$1@fred.mathworks.com>,
> evan klinger <eklinger1-nospam@cox.net> wrote:
> >roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
> >message <fhaiam$k8$1@canopus.cc.umanitoba.ca>...
> >> In article <fhag57$bca$1@fred.mathworks.com>,
> >> evan klinger <eklinger1-nospam@cox.net> wrote:
> >> >Is there a way to increase the 25,000 character limit? I
> >> >really need to solve up to cv12,
>
> >> you are trying to find the closed form solution to a
> >polynomial
> >> of order 12 that does not happen to be factorable down
>
> >So what are my options here? I would like a closed form
> >solution so that I can plug in values into the variables
> >directly and solve, but if that's not possible, is there any
> >alternative?
>
> What kind of ranges of values can your constants take?
> I've been playing around with some test values, and it
looks to
> me as if you could easily encounter severe problems with
loss of
> precision in your calculations, especially if cr > 2.
>
> With random values between -1 and 1 for the cv12 variables, it
> was fairly common in my trials for there to be an x root
somewhere
> just a bit more than 2, but that by x = 3 that the
function value
> had risen beyond 10^8 or 10^9 -- but with another one of
my trials,
> the roots ranged from -7 to +4 .
>
> Unless your values (especially cr) are fairly small, I suspect
> you will lose precision so quickly as to make it very
doubtful that
> you would be able to accurately find the equation roots.
>
> cr is the biggest contribution to the problem: as well as
> cv12 being a 12th order polynomial in x, it is also a 12th
order
> polynomial in cr (cr and x will have the same polynomial order
> as each other in each cv*). Thus even a small change in cr can
> make a big change in the x roots. And a small change in cr can
> result in -different- 12th order roots becoming the
real-valued
> roots, so even if you had an closed form solution for one
of the
> 12th-order roots, that root might become imaginary with a
small
> change in parameters.
>
> For example, in a particular trial I did
> with random values for the parameters except p3, for p3=1/3,
> the first and 7th roots were the real roots, but for p3=1/2,
> the 6th and 7th roots were the real roots -- and p3
appears only
> in its first power in cv12.
>
>
> If I recall correctly, I have seen mention of a FEX
contribution
> that generalizes fzero() to find multiple roots. You are
> probably going to have to use something like that, and to
> figure out what you want to do with the multiple roots
that are
> found, but with some kinds of checks in place to ensure
that you
> are not finding spurious roots (or failing to find valid
roots)
> due to loss of precision.
> --
> We regret to announce that sub-millibarn resolution
bio-hyperdimensional
> plasmatic space polyimaging has been delayed until the release
> of Windows Vista SP2.

cr will always be between 0 and 1
db will be in the millions
p is usually between 0 and 1,000,000
cv[n] is usually between 0 and a few million

This is related to cost of insurance for a life insurance policy

In article <fhd1pu$j8d$1@fred.mathworks.com>,
evan klinger <eklinger1-nospam@cox.net> wrote:
>roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
>message <fhb2v2$jqi$1@canopus.cc.umanitoba.ca>...

>> What kind of ranges of values can your constants take?
>> I've been playing around with some test values, and it
>looks to
>> me as if you could easily encounter severe problems with
>loss of
>> precision in your calculations, especially if cr > 2.

>cr will always be between 0 and 1
>db will be in the millions
>p is usually between 0 and 1,000,000
>cv[n] is usually between 0 and a few million

cv12 has no real solutions if cv0 > db, which is a possibility
under the constraints you have given (db "in the millions"
could be as low as 1 million; cv0 up to "a few million"
could exceed that value.)

I tried each combination of the lower and upper value you gave
for each parameter. Taking only the extremes, there were
real roots for each case unless cv0 > db. The cases that had
only one real root, that root was >= 0. The cases that had
two real roots, the second root was < 0.

In random trials over the range of parameters, -usually- there was
one positive and one negative root, but in one of the cases there
were only two negative roots. In one of the cases, the positive root
was about +46 instead of the typical 0.08 - 0.25 range. That's going
to make it difficult to give good starting guesses for the ranges,
since your precision will usually be completely washed out before x = 4.
--
We regret to announce that sub-millibarn resolution bio-hyperdimensional
plasmatic space polyimaging has been delayed until the release
of Windows Vista SP2.

roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
message <fhe9fk$l62$1@canopus.cc.umanitoba.ca>...
> In article <fhd1pu$j8d$1@fred.mathworks.com>,
> evan klinger <eklinger1-nospam@cox.net> wrote:
> >roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
> >message <fhb2v2$jqi$1@canopus.cc.umanitoba.ca>...
>
> >> What kind of ranges of values can your constants take?
> >> I've been playing around with some test values, and it
> >looks to
> >> me as if you could easily encounter severe problems with
> >loss of
> >> precision in your calculations, especially if cr > 2.
>
> >cr will always be between 0 and 1
> >db will be in the millions
> >p is usually between 0 and 1,000,000
> >cv[n] is usually between 0 and a few million
>
> cv12 has no real solutions if cv0 > db, which is a possibility
> under the constraints you have given (db "in the millions"
> could be as low as 1 million; cv0 up to "a few million"
> could exceed that value.)
>
> I tried each combination of the lower and upper value you gave
> for each parameter. Taking only the extremes, there were
> real roots for each case unless cv0 > db. The cases that had
> only one real root, that root was >= 0. The cases that had
> two real roots, the second root was < 0.
>
> In random trials over the range of parameters, -usually-
there was
> one positive and one negative root, but in one of the
cases there
> were only two negative roots. In one of the cases, the
positive root
> was about +46 instead of the typical 0.08 - 0.25 range.
That's going
> to make it difficult to give good starting guesses for the
ranges,
> since your precision will usually be completely washed out
before x = 4.
> --
> We regret to announce that sub-millibarn resolution
bio-hyperdimensional
> plasmatic space polyimaging has been delayed until the release
> of Windows Vista SP2.

This is what I'm currently using to solve it. It takes about
10 seconds to solve. Is there a way to optimize it?

Thanks

function [coi] = calc_coi(cv0, db, cr,
p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,cvgoal)

syms x
cv1 = (cv0+p1-x*(db-cv0))*(1+cr);
simplify(cv1);
cv2 = (cv1+p2-x*(db-cv1))*(1+cr);
simplify(cv2);
cv3 = (cv2+p3-x*(db-cv2))*(1+cr);
simplify(cv3);
cv4 = (cv3+p4-x*(db-cv3))*(1+cr);
simplify(cv4);
cv5 = (cv4+p5-x*(db-cv4))*(1+cr);
simplify(cv5);
cv6 = (cv5+p6-x*(db-cv5))*(1+cr);
simplify(cv6);
cv7 = (cv6+p7-x*(db-cv6))*(1+cr);
simplify(cv7);
cv8 = (cv7+p8-x*(db-cv7))*(1+cr);
simplify(cv8);
cv9 = (cv8+p9-x*(db-cv8))*(1+cr);
simplify(cv9);
cv10 = (cv9+p10-x*(db-cv9))*(1+cr);
simplify(cv10);
cv11 = (cv10+p11-x*(db-cv10))*(1+cr);
simplify(cv11);
cv12 = (cv11 + p12 - (x)*(db-cv11))*(1+cr);
simplify(cv12);
eq = cv12 - cvgoal;
ret = solve(eq);

coi = ret(1);

If I remove the simplify() commands it takes even longer, up
to 18 seconds.

"evan klinger" <eklinger1-nospam@cox.net> wrote in message
<fhf89b$srr$1@fred.mathworks.com>...
> roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
> message <fhe9fk$l62$1@canopus.cc.umanitoba.ca>...
> > In article <fhd1pu$j8d$1@fred.mathworks.com>,
> > evan klinger <eklinger1-nospam@cox.net> wrote:
> > >roberson@ibd.nrc-cnrc.gc.ca (Walter Roberson) wrote in
> > >message <fhb2v2$jqi$1@canopus.cc.umanitoba.ca>...
> >
> > >> What kind of ranges of values can your constants take?
> > >> I've been playing around with some test values, and it
> > >looks to
> > >> me as if you could easily encounter severe problems with
> > >loss of
> > >> precision in your calculations, especially if cr > 2.
> >
> > >cr will always be between 0 and 1
> > >db will be in the millions
> > >p is usually between 0 and 1,000,000
> > >cv[n] is usually between 0 and a few million
> >
> > cv12 has no real solutions if cv0 > db, which is a
possibility
> > under the constraints you have given (db "in the millions"
> > could be as low as 1 million; cv0 up to "a few million"
> > could exceed that value.)
> >
> > I tried each combination of the lower and upper value
you gave
> > for each parameter. Taking only the extremes, there were
> > real roots for each case unless cv0 > db. The cases that had
> > only one real root, that root was >= 0. The cases that had
> > two real roots, the second root was < 0.
> >
> > In random trials over the range of parameters, -usually-
> there was
> > one positive and one negative root, but in one of the
> cases there
> > were only two negative roots. In one of the cases, the
> positive root
> > was about +46 instead of the typical 0.08 - 0.25 range.
> That's going
> > to make it difficult to give good starting guesses for the
> ranges,
> > since your precision will usually be completely washed out
> before x = 4.
> > --
> > We regret to announce that sub-millibarn resolution
> bio-hyperdimensional
> > plasmatic space polyimaging has been delayed until the
release
> > of Windows Vista SP2.
>
> This is what I'm currently using to solve it. It takes about
> 10 seconds to solve. Is there a way to optimize it?
>
> Thanks
>
> function [coi] = calc_coi(cv0, db, cr,
> p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,cvgoal)
>
> syms x
> cv1 = (cv0+p1-x*(db-cv0))*(1+cr);
> simplify(cv1);
> cv2 = (cv1+p2-x*(db-cv1))*(1+cr);
> simplify(cv2);
> cv3 = (cv2+p3-x*(db-cv2))*(1+cr);
> simplify(cv3);
> cv4 = (cv3+p4-x*(db-cv3))*(1+cr);
> simplify(cv4);
> cv5 = (cv4+p5-x*(db-cv4))*(1+cr);
> simplify(cv5);
> cv6 = (cv5+p6-x*(db-cv5))*(1+cr);
> simplify(cv6);
> cv7 = (cv6+p7-x*(db-cv6))*(1+cr);
> simplify(cv7);
> cv8 = (cv7+p8-x*(db-cv7))*(1+cr);
> simplify(cv8);
> cv9 = (cv8+p9-x*(db-cv8))*(1+cr);
> simplify(cv9);
> cv10 = (cv9+p10-x*(db-cv9))*(1+cr);
> simplify(cv10);
> cv11 = (cv10+p11-x*(db-cv10))*(1+cr);
> simplify(cv11);
> cv12 = (cv11 + p12 - (x)*(db-cv11))*(1+cr);
> simplify(cv12);
> eq = cv12 - cvgoal;
> ret = solve(eq);
>
> coi = ret(1);
>
> If I remove the simplify() commands it takes even longer, up
> to 18 seconds.
>


Bump