mm-929
===
Subject: Solve[] and Eliminate[] choke on simple systems of equations
I am trying to convert some of my old stuff from another system to
Mathematica.
Below is an example of a Markov model that takes the other system a couple
of seconds to solve, but Mathematica chews on forever (I gave up after
about 30 minutes).
Since I'm a Mathematica newbie, perhaps there is some limitation I am
not aware of?
eqns = {3*L*P333 == U*(P323 + P322 + P223 + P221 + P120),
(K + 2*L + U)*P323 == 3*L*P333,
(K + 2*U + L)*P313 == 2*L*P323,
(2*L + U)*P322 == 2*U*P312,
(2*U + L)*P312 == 2*L*P322 + 3*U*P302,
3*U*P302 == L*P312 + L*P313,
(2*L + U)*P223 == U*(2*P213 + P212) + 2*U*P313 + U*(2*P113 + P111) +
K*P323,
(K + 2*U + L)*P213 == 2*L*P223,
(2*U + L)*P212 == L*P221 + 2*U*P202,
3*U*P202 == L*(P213 + P212 + P211),
(2*L + U)*P221 == U*(P212 + 2*P211),
(2*U + L)*P211 == L*P221 + U*P202
(L + 2*U)*P113 == K*P313 + K*P213 + U*P102,
3*U*P102 == L*P113 + L*P111,
(2*U + L)*P111 == 2*U*P102 + 2*L*P120,
(2*L + U)*P120 == U*P111,
P333 + P323 + P313 + P322 +
P312 + P302 + P223 + P213 + P212 + P202 + P221 + P211 + P113 +
P102 +
P111 + P120 == 1}
vars = {P333, P323, P313, P322, P312, P302, P223, P213, P212, P202,
P221,
P211, P113, P102, P111, P120}
A = Solve[eqns, vars]
===
Subject: Re: Solve[] and Eliminate[] choke on simple systems of equations
you have 17 equations and 16 variables! I think there is something fishy
with the tripple equation:
(2*U + L)*P211 == L*P221 + U*P202 (L + 2*U)*P113 == K*P313 + K*P213 +
U*P102
have fun Daniel
===
Subject: Re: Integrate (a curious result)
I believe that the real reason is very close to what you have guessed
on your post!
More or less it must be something on the same vein.
Otherwise it could be an internal problem of Integration
agorithm for Mathematica 5.2 for Windows (the integral is
platform dependent).
Anyway...
Dimitris
=CF/=C7 Michael Weyrauch =DD=E3=F1=E1=F8=E5:
ite integration. So, please, keep on asking these questions...
see what is going on,
knows the source
lgorithms for
discontinuities along
hich may go wrong
an idea how I should
way as you do in your post.
if the various limits exist or not.
ess that it is this what happens in
out the limits at this singularity and therefore
onsistent with your observation
culates the limits at the endpoints of the
tegral is Infinity.
it actually works...
===
Subject: Re: Efficient BesselJ and LegendreP Evaluation
Hi Antonio.
Your post apperared in quite unreadable format (take a look at the
functions arguments).
Anyway from our private communication I was able to see your integrals
you are interested in.
I will what I can do.
Dimitris
=CF/=C7 Antonio =DD=E3=F1=E1=F8=E5:
Try
^4);
);
E=B1Lim = =CE=B1Ma=
er))^2];
CE=B1] +
CE=B1] +
Exp[-(f*(S=
E=B1]]*(Bess=
=B1])) -
+
B1]]*LegendreP=
Legendre[n, m,
=CF=81o*Si=
===
Subject: Re: Efficient BesselJ and LegendreP Evaluation
Sorry for the greek letters in the code below I have replaced them
with their names.
Clear[Global`*];
Off[General::spell];
Off[General::spell1];
lambda = 0.8;
w = 2500;
NA = 1.25;
f = 1700;
d = 50.;
nGlass = 1.4883 + 1.9023/(10^2*lambda^2) - 4.25016/(10^3*lambda^4);
nWater = 1.3253 + 2.7553/(10^3*lambda^2) + 3.779/(10^5*lambda^4);
nAir = 1.0001;
alphaMax = ArcSin[NA/nWater];
ko = (2*Pi*nAir)/lambda;
Nmax = 20;
GridNM = Table[{n, m}, {n, 1, Nmax}, {m, 0, n}];
alphaCrit = ArcSin[nWater/nGlass];
If[alphaCrit < alphaMax, alphaLim = alphaCrit, alphaLim = alphaMax];
cosA2[(alpha_)?NumericQ] := Sqrt[1 - (nGlass*(Sin[alpha]/nWater))^2];
ts[(alpha_)?NumericQ] := 2*nGlass*(Cos[alpha]/(nGlass*Cos[alpha] +
nWater*cosA2[alpha]));
tp[(alpha_)?NumericQ] := 2*nGlass*(Cos[alpha]/(nWater*Cos[alpha] +
nGlass*cosA2[alpha]));
DBessel[m_, x_] := (1/2)*(BesselJ[-1 + m, x] - BesselJ[1 + m, x]);
DLegendre[n_, m_, x_] := ((-1 - n)*x*LegendreP[n, m, x] + (1 - m +
n)*LegendreP[1 + n, m, x])/(-1 + x^2);
ITM[(n_Integer)?Positive, (m_Integer)?NonNegative, zo_Real, rhoo_Real] :=
Block[{func}, func[(alpha_)?NumericQ] :=
Sqrt[Cos[alpha]]*Exp[-(f*(Sin[alpha]/w))^2]*Exp[(-I)*ko*nWater*zo*cosA2[alph
a]]*
Exp[(-I)*ko*d*(nGlass*Cos[alpha] -
nWater*cosA2[alpha])]*(m^2*ts[alpha]*LegendreP[n, m,
cosA2[alpha]]*(BesselJ[m,
ko*nGlass*rhoo*Sin[alpha]]/(ko*nGlass*rhoo*Sin[alpha])) -
tp[alpha]*(nGlass/nWater)^2*Sin[alpha]^2*DBessel[m,
ko*nGlass*rhoo*Sin[alpha]]*DLegendre[n, m, cosA2[alpha]] +
I*m*(ts[alpha]*DBessel[m,
ko*nGlass*rhoo*Sin[alpha]]*LegendreP[n, m, cosA2[alpha]] -
tp[alpha]*(nGlass/nWater)^2*Sin[alpha]^2*DLegendre[n, m,
cosA2[alpha]]*
(BesselJ[m,
ko*nGlass*rhoo*Sin[alpha]]/(ko*nGlass*rhoo*Sin[alpha]))));
1000]];
zo = Random[Real, {-10., 10.}];
rhoo = Random[Real, {0., 10.}];
Timing[MI1 = Apply[ITM[#1, #2, zo, rhoo] & , GridNM, {2}];]
===
Subject: Timing in Mathematica
I am using Timing[.......] to calculate the time taken by an
expression in Mathematica. But it is not showing values below 0.01
seconds. For values below 0.01s it just shows 0. Second
Is it possible to get the value of time lower than 0.01s by somehow
lowering the least count. I am using Mathematica 5.2 on Windows. I was
getting lower values of time on Mathematica(5.2) on a Linux computer.
Amit
===
Subject: Re: Timing in Mathematica
It is very easy.
Let for example
foo = {x^n, Cos[x], Log[x]/x, Log[Sin[x]^2]*Tan[x], BesselJ[n, x],
(2*x)/((x + 1)*(x^3 + 3*x^2 + 2*x + 1))}
{x^n, Cos[x], Log[x]/x, Log[Sin[x]^2]*Tan[x], BesselJ[n, x], (2*x)/((1
+ x)*(1 + 2*x + 3*x^2 + x^3))}
Then
InputForm[(Timing[Integrate[#1, x]] & ) /@ foo]
{{0.031*Second, x^(1 + n)/(1 + n)}, {0.015999999999999986*Second,
Sin[x]}, {0.*Second, Log[x]^2/2},
{0.9530000000000001*Second, Log[Sec[x/2]^2]^2 + 2*Log[Sec[x/
2]^2]*Log[(Cos[x]*Sec[x/2]^2)/2] +
Log[Sec[x/2]^2]*Log[Sin[x]^2] - 2*Log[Sec[x/2]^2]*Log[-1 + Tan[x/
2]^2] - Log[Sin[x]^2]*Log[-1 + Tan[x/2]^2] +
Log[Tan[x/2]^2]*Log[-1 + Tan[x/2]^2] + 2*PolyLog[2, Sec[x/2]^2/2]
+ PolyLog[2, Cos[x]*Sec[x/2]^2] +
PolyLog[2, -Tan[x/2]^2]}, {0.030999999999999917*Second, 2^(-1 -
n)*x^(1 + n)*Gamma[(1 + n)/2]*
HypergeometricPFQRegularized[{(1 + n)/2}, {1 + n, (3 + n)/2}, -
(x^2/4)]},
{0.016000000000000014*Second, 2*(-Log[1 + x] + RootSum[1 + 2*#1 +
3*#1^2 + #1^3 & ,
(Log[x - #1] + 2*Log[x - #1]*#1 + Log[x - #1]*#1^2)/(2 + 6*#1 +
3*#1^2) & ])}}
Dimitris
=CF/=C7 amitsoni.1984@gmail.com =DD=E3=F1=E1=F8=E5:
===
Subject: Re: Timing in Mathematica
Hi Amit,
if it is important to get these timings (because the functions are called
in
loops), I would do it with
and divide the result by 10^4.

===
Subject: Re: Timing in Mathematica
amitsoni.1984@gmail.com =EDrta:
After doing some experimentation, I found that on my system
AbsoluteTiming returns integer multiples of 0.015625 seconds. Timing
returns the same values, but rounded to $TimeValue of 0.001.
(Mathematica 5.2 on WinXP SP2)
But you still shouldn't assume that timings are precise to 0.015625.
For some reason, when I record several timing values with Table, the
first one is always larger then the rest (significantly larger for
short times).
Here are two examples with very short timings:
Table[First@
AbsoluteTiming[Do[NSum[i,{i,1,1000}],{1000}]]/Second,{100}]/
0=2E015625
{36.,25.,26.,27.,29.,29.,28.,28.,27.,27.,26.,26.,26.,26.,25.,26.,25.,26.,26=

.,
25.,26.,25.,26.,25.,26.,28.,25.,27.,26.,27.,26.,25.,27.,26.,26.,25.,26.,27.=

,
26.,26.,25.,26.,27.,25.,26.,25.,25.,28.,26.,26.,28.,25.,26.,25.,27.,26.,26.=

,
26.,26.,26.,25.,25.,26.,26.,29.,25.,26.,25.,25.,26.,25.,26.,25.,25.,26.,27.=

,
25.,26.,25.,26.,25.,25.,26.,25.,25.,26.,25.,25.,26.,25.,25.,26.,25.,26.,25.=

,
25.,25.,26.,25.,26.}
Table[First@
AbsoluteTiming[Do[NSum[i,{i,1,1000}],{100}]]/Second,{100}]/
0=2E015625
{6.,6.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.=

,3.
,
3=2E,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,3.,2=

.,3.
,
2=2E,3.,2.,3.,2.,3.,2.,3.,2.,3.,3.,2.,3.,2.,3.,2.,3.,3.,2.,3.,2.,3.,2.,3.,2=

.,3.
,3.,2.,3.,2.,3.,2.,3.,2.,3.,2.,3.,3.,2.,3.,2.,3.,3.,2.,3.,2.,3.,2.}
Other calculations inside Do[] or using Timing instead of
AbsoluteTiming give similar results.
I hope that this will persuade you not try to measure very short
execution times.
===
Subject: Visual Basic
I'm a french student and I use mathematica, I just want to know if it is
po=
ssible to run Mathematica from a Visual Basic program and to return the
mat=
hematica output in my window application.
_________________________________________________________________
Essayez Live.com et cr=E9ez l'Internet qui vous ressemble : infos, sports,
=
m=E9t=E9o et bien plus encore !
http://www.live.com/getstarted
===
Subject: Re: Depurating tool
I think you need Wolfram Workbench, a great new tool. Unfortunately, I
guess,
it is only available right now to Premier Service customers. But that you
need
to ask customer support of WRI.
Michael Weyrauch
===
Subject: Re: verification
...
Hi Dimitris,
try this one:
Reduce[And @@ Equal @@@ Subsets[foo, {2}] // ComplexExpand]
Peter
===
Subject: Re: verification
OK, so you easily showed that foo[[2]] and foo[[3]] are the same. Here's
the easiest way I know to show that foo[[1]] and foo[[3]] are the same:
In[4]:= FullSimplify[TrigToExp[foo[[1]]]] == FullSimplify[foo[[3]]]
Out[4]= True
What really disturbs me is that I can also show that foo[[1]] and
foo[[3]] are NOT the same:
In[6]:= FullSimplify[TrigToExp[foo[[1]]] == foo[[3]]]
messages regarding Internal precision limit snipped
Out[6]= False
It seems that Mathematica is asserting that two _equal_ expressions
are _not equal_. What am I missing?! (Of course, if Mathematica had merely
left the logical expression unevaluated, I wouldn't have been disturbed...)
BTW, I haven't found a direct way to use Mathematica to show that foo[[1]]
and foo[[2]] are the same.
David W. Cantrell
===
Subject: Re: verification
In fact it is quite easy to prove using Mathematica that all the
elements of foo must be equal although it needs a little human
intervention.
foo = {ArcTan[8/(1 - Sqrt[-15 - 4*I])] + ArcTan[8/(1 + Sqrt[-15 -
4*I])] + ArcTan[8/(1 - Sqrt[-15 + 4*I])] +
ArcTan[8/(1 + Sqrt[-15 + 4*I])], ArcTan[3] + ArcTan[5] +
ArcTan[41/3] + ArcTan[21], 2*Pi - ArcTan[1/4] - ArcTan[5/12]};
FullSimplify[Tan[foo[[2]]]-Tan[foo[[3]]]]
0
Since both foo[[2]] and foo[[3]] are real, in order to have equal
values of Tan they must differ by a multiple of Pi but:
Abs[foo[[2]]-foo[[3]]]<Pi
True
Essentially in the same way we can prove that foo[[1]] == foo[[2]].
Mathematica does not seem to be able to verify symbolically that foo
[[1]] is real, however, from defintion of Tan we know that it is
impossible for a non-real and real numbers to have equal values of
Tan, but
FullSimplify[Tan[foo[[1]]]-Tan[foo[[2]]]]
0
so foo[[1]] is also real. Exactly as above, we only need to check that
Abs[foo[[1]]-foo[[2]]]<Pi
True
That shows that they are the elements of the list foo are equal real
numbers.
Andrzej Kozlowski
===
Subject: Re: Kernel closes after nested function returns
Mathematica copy the data quite often and I would expect, that you
run out of memory or the kernel need more memory than the operating
system allows. You may hinder the copy process by using a Block[] in
func A and use in func B the variable name of func A
funcA[listInput_]:=Block[{aVeryLargeList=listInput},
...
funcB[]
...
]
funcB[]:=Module[{someVariable},
...
aVeryLargeList*someVariable
...
]
eventual with the full context.
Or you may transfer the symbol with the HoldFirst attribute
or you may use a 64 bit system.
Jens
===
Subject: Re: Re: Cantor Function problem
The corrected code I gave you did work on my system and it will works
on yours if you enter it (or cut and paste it) *verbatim*. In other
words, you must respect the case of the characters for Mathematica
distinguishes a c (lowercase c) from a C (uppercase c).
?Connect
link listening on a named port.
Connect is a system function and Mathematica complains about your
attempt to redefine it. Use connect (lowercase only) as given in my
previous message.
It might be worse reading, or at least browsing, the first part of
_The Mathematica Book_ to avoid such pitfails [1].
HTH,
Jean-Marc
[1] Part 1: A Practical Introduction to Mathematica, _The Mathematica
Book_,
http://documents.wolfram.com/mathematica/book/section-1
===
Subject: Re: Cantor Function problem
I apologize for this nonsense. What I meant to say was just the
opposite! The correct sentence reads It might be worth reading, or at
least browsing, the first part of _The Mathematica Book_ to avoid such
pitfalls.
Jean-Marc
===
Subject: Re: Urgent( Problem with Mathematica-5)
The corrected code I gave you did work on my system and it will works
on yours if you enter it (or cut and paste it) *verbatim*. In other
words, you must respect the case of the characters for Mathematica
distinguishes a c (lowercase c) from a C (uppercase c).
?Connect
link listening on a named port.
Connect is a system function and Mathematica complains about your
attempt to redefine it. Use connect (lowercase only) as given in my
previous message.
It might be worse reading, or at least browsing, the first part of
_The Mathematica Book_ to avoid such pitfalls [1].
HTH,
Jean-Marc
[1] Part 1: A Practical Introduction to Mathematica, _The Mathematica
Book_,
http://documents.wolfram.com/mathematica/book/section-1
===
Subject: Re: Re: comments in import files
Hi Dave,
a set of tools. Checking their web site, I found plenty of information
about Windows Services for Unix (SFU) 3.5, which is too fatty to my
taste :-) I will stick with Cygwin on Windows and otherwise genuine
Unix and Linux systems.
Jean-Marc
===
Subject: Re: Simplification with Integers assumption
I tried in Math 4 (WinXP) your expression:
In[1]:= FullSimplify[ Sum[ n! / (2*p)! / (n - 2*p)! , {p, 1, Infinity}]

Out[1]= -1 + 2^(-1+n)
In[2]:= Simplify[ Sum[ n! / (2*p)! / (n - 2*p)! , {p, 1, Infinity}] , {n

n [Element] Integers}]
Out[2]= (-2+2^n) n! / (2 Gamma[1+n])
Math 5.2 makes me puzzled!
did
===
Subject: Re: Simplification with Integers assumption
The problem is that
Assuming[Element[n,Integers],Simplify[Sin[n Pi]]]
0
But
Simplify[Gamma[n], n =E2=88=88 Integers]
Gamma[n]
(even adding the assumption n<0 will still return Gamma[n]).
Since
Sum[n!/(2*p)!/(n - 2*p)!, {p, 1, Infinity}]
-(((-2 + 2^n)*n!*Gamma[-n]*Sin[n*Pi])/(2*Pi))
simplifying with the assumption ELement[n,Integers] will result in
Gamma[-n]*0 and since in Mathematica any symbol s (except explicit
Infitiniy and its variants) s*0 is always 0, you will get 0.
The best way to avoid the problem is to try to get an expression for
the sum not involving the Gamma function at all, if possible. In this =
since FullSimplify without the condition Element[n,Integers] returns
an expression not involving Gamma, you cant then Simplify again using =
the this assumption, although in this particular case it won't make
any difference. This is one of the cases where one needs to pay
careful attention to the effect of different assumptions in
FullSimplify may have on the the expression they are being applied
to. The danger is that if one of the assumptions implies that one of
the factors in in a product is zero the whole product will be reduced =
to 0, unless one of the other factors are explicitly known to
Mathematica=E2=80=99s evaluator be infinite under the given assumptions. =
An
example of the latter case is:
FullSimplify[Sum[n!/(2*p)!/(n - 2*p)!, {p, 1, Infinity}], {n == 3}]
While this partiular problem could perhaps be solved by adding the
knowledge that Gamma[n] is CompexInfinity for any negative integer to =
Simplify and FullSimplify some problems of this kind will always
remain since Mathematica's evaluator turns x*0 to 0 even if nothing
is known about x (and in particular when x is an implicit but not
explicit Infinity).
Andrzej Kozlowski
===
Subject: Re: Simplification with Integers assumption
works o.k. in my version:
0] and Simplify[Sum[n!/(2*p)!/(n - 2*p)!, {p, 1, Infinity}], Assumptions
-(((-2 + 2^n)*n!*Gamma[-n]*Sin[n*Pi])/(2*Pi))
and with FullSimplify both give the simpler:-1 + 2^(-1 + n)
Daniel
===
Subject: Re: Simplification with Integers assumption
........................Gamma[-n] is what you intended above, I believe.
Well, it certainly does not surprise me.
But your post did lead me to find something surprising. Suppose we fix n as
a positive integer -- say, 3 for example:
In[6]:= Sum[3/((2p)!(3 - 2p)!), {p, 1, Infinity}]
results in nothing but warning messages and gibberish. OTOH, if we replace
Infinity with a large integer -- say, 100 or 1000 -- Mathematica has no
trouble giving the correct answer.
In[7]:= Sum[3/((2p)!(3 - 2p)!), {p, 1, 1000}]
Out[7]= 3/2
Of course, that answer is _also_ the correct one for the infinite sum at
which Mathematica balked. Why couldn't Mathematica do the infinite sum? Is
this problem fixed in the development version perhaps?
David W. Cantrell
===
Subject: Series and order
I wish to expand the expression e1, below, in a Taylor's series
expansion about z0 in the two variables z1 and z2 and include only
first order terms (c1, c2, a and b are coefficients). Terms in z1 and
z2 are of similar order so that products of their first order terms
are second order and should be discounted compared to first order
terms. Problem 1: how do I tell Series to do this? Problem 2: the
coefficients c1 and c2 are of similar order to z1-z0 and thus
additional terms can be set to zero as second order. How do I do
this?
Below I give a warm up problem for a function f[z1,z2] which just
seems to expand in terms of one variable and then in terms of the
other without regard for term order. If we can solve the first problem
then it may be possible to solve the second by expanding in terms of
the coefficients as well.
Hugh Goyder
Series[f[z1, z2], {z2, z0, 1}, {z1, z0, 1}]
SeriesData[z2, z0, {SeriesData[z1, z0,
{f[z0, z0], Derivative[1, 0][f][z0, z0]}, 0, 2, 1],
SeriesData[z1, z0, {Derivative[0, 1][f][z0, z0],
Derivative[1, 1][f][z0, z0]}, 0, 2, 1]}, 0, 2, 1]
e1 = (z0^2 - 2*I*c1*z0*z1 - z1^2)*(z0^2 - 2*I*c2*z0*z2 -
z2^2) - 4*z1*z2*(I*a*z0 + b*Sqrt[z1*z2])^2;
Simplify[Normal[Series[e1, {z2, z0, 1}, {z1, z0, 1}]],
{0 < z0}]
2*z0^2*((2 - I*a*b - 2*b^2)*z0^2 +
(2*a^2 - 9*I*a*b - 8*b^2 - 2*(-I + c1)*(-I + c2))*z1*z2 +
z0*((-2 + 3*I*a*b + 4*b^2 - 2*I*c1)*z1 +
(-2 + 3*I*a*b + 4*b^2 - 2*I*c2)*z2))
===
Subject: Re: choke on simple systems
hurt (it was a direct reading of the Markov model).
----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.
===
Subject: ftp from mathematica code
is it possible to ftp (get) a file from within either/both a
mathematica notebook code or/and an m file code running on the command
line (should be equivalent, right?).
thx for any help.
===
Subject: summation is equivalent to dot product??
Comparing the caluclations done by summation and dot product of the
matrix, it turns out that they are different, which i initially
thought to be equal. what is the
equivalent summation form for the dot product??
In[29]:=
F[p_,q_,r_,s_]:=p+2q+r-s;
In[31]:=
MatF=Table[F[p,q,r,s],{p,4},{q,4},{r,4},{s,4}];
In[32]:=
c={1,2.1,3,4};
In[37]:=
c=2EMatF.c//MatrixForm
sumf[
p_, q_] := =E2=88=91r=E2=88=91s F[p, q, r, s]*c[[r]] c[[s]];
In[41]:=
Table[sumf[p,q],{p,4},{q,4}]//MatrixForm
===
Subject: Re: Finding unknown parameters using Mathematica
I've been looking the system of 5 equations you quoted. I don=B4t know
if it has a solution. Numerically ,I guess the function that could
work in your case is FindRoot. You have to give initial values to the
variables, e.g. all ones. I think FindRoot minimizes the absolute
value of the l.h.s. of your equations and determines an approximate
solution. The problem is that the first two equations have r1 or r2 in
the denominators, so, the numerical search tends to increase these two
values without limit. Then I multiplied the first equation by r1 and
the 2nd by r2 and using FindMinimum in order to minimize the sums of
the squares of the l.h.s.'s of your equations (increasing the number
of MaxIterations) got this:
r1 0.000090274767504370
r2 5.520801180767231e - 6
b1 2.4543955633519086
b2 3.7417707839019467
p 864.1266904291423
That makes the first four equations approx. 10e-10 and the fifth 10e-7.
The system is very non-linear though....
I don't know. I hope this could help.
I also suggest not to put directly the Mathematica code in the
message because there are often problems interpreting it.
===
Subject: Finding unknown parameters using Mathematica
I'm very new to mathematica and trying to solve a set of non-linear system
of equations to find the unknown parameters for a bivariate distribution. =
I've 5 unknown parameters (i.e. b1,b2,r1,r2,p) and 5 set of equations. I =
tried to get the general solution , but could not. I don't know how to =
solve these nonlinear equations to get the unknown parameters. And also =
not aware how to give initial value in the solve function or any other =
function (e.g all parameters b1,b2,r1,r2,p if I give initial boundary =
value =1).
In the below codes, x1 and x2 are 2 variables (e.g. data from 2 quality =
characteristics; stress and strain etc. etc.)
Any comments / suggestions how to solve these equation for b1,b2,r1,r2,p =
where as putting intial value for all these unknown parameters =1
Ahmad S.
n=4
x1={1,2,3,4}
x2={1.7,3.8,4.9,4.6}
4
{1,2,3,4}
{1.7,3.8,4.9,4.6}
!(*
RowBox[{
RowBox[{
StyleBox[L,
StyleBox[=,
StyleBox[(n*Log[p] + n*Log[p + 1] + n*
Log[b1] + n*Log[r1] + n*Log[b2] + n*Log[r2] + ((b1 - 1)) =
(=E2=88=91+(j =
1)%n Log[x1[([
j])]]) + ((b2 -
1)) (=E2=88=91+(j = 1)%n Log[
x2[([j])]]) - ((p + 2)) (=E2=88=91+(j = 1)%n
Log[1 +
r1*((x1[([j])]^b1)) +
r2*((x2[([j])]^b2))])),
StyleBox[ ,
!(4.980920826406141` (((-1) + b2)) + (((-1) + b1))
((
Log[2] + Log[3] + Log[
4])) + 4 Log[b1] + 4 Log[b2] + 4 Log[
p] + 4 Log[1 + p] + 4 Log[r1] + 4 Log[r2] - ((2 + p))
((
Log[1 + r1 + 1.7`^b2
r2] + Log[1 + 2^b1 r1 + 3.8`^b2 r2] +
Log[1 + 4^b1 r1 + 4.6`^b2 r2] + Log[1 + 3^b1 r1 + =
4.9`^b2
r2])))
Eqn1=D[L,r1][Equal]0
!(4/r1 - ((2 +
p)) ((1/(1 + r1 + 1.7`^b2
r2) + 2^b1/(1 + 2^b1 r1 + 3.8`^b2 r2) +
4^b1/(1 +
4^b1 r1 + 4.6`^b2 r2) + 3^b1/(1 + 3^b1 r1 +
=
4.9`^b2
r2))) [Equal] 0)
Eqn2=D[L,r2][Equal]0
!(4/r2 - ((2 +
p)) ((1.7`^b2/(1 +
r1 + 1.7`^b2 r2) + 3.8`^b2/(1 + 2^b1 r1 +
3.8`^b2 r2) =
+
4.6`^b2/(1 +
4^b1 r1 + 4.6`^b2 r2) + 4.9`^b2/(1 + 3^b1 r1
+ =
4.9`^b2
r2))) [Equal] 0)
Eqn3=D[L,b1][Equal]0
!(4/b1 + Log[2] + Log[3] + Log[4] - ((
2 + p)) (((2^b1 r1
Log[2])/(1 + 2^b1 r1 + 3.8`^b2
r2) + (3^b1 r1 Log[3])/(1 + 3^b1 r1 +
4.9`^b2 r2) + =

(4^b1 r1 Log[4])/(1 + 4^b1 r1 + 4.6`^b2 r2)))
[Equal] 0)
Eqn4=D[L,b2][Equal]0
!(((4.980920826406141`)([InvisibleSpace])) +
4/b2 - ((2 + p)) (((0.5306282510621704` 1.7`^b2
r2)/(1 =
+
r1 + 1.7`^b2
r2) + (1.33500106673234` 3.8`^b2 r2)/(1 +
2^b1 r1 + =

3.8`^b2 r2) + (1.5260563034950492` 4.6`^b2 r2)/(1 +
4^b1 r1 + =

4.6`^b2 r2) + (1.589235205116581` 4.9`^b2 r2)/(1 +
3^b1 r1 + =
4.9`
^b2 r2))) [Equal] 0)
Eqn5=D[L,p][Equal]0
!(4/p + 4/(1 + p) - Log[1 + r1 + 1.7`^b2
r2] - Log[1 + 2^b1 r1 + 3.8`^b2 r2] - Log[1 + 4^b1
r1 + 4.6`^b2 r2] - Log[1 + 3^b1 r1 + 4.9`^b2 r2] [Equal]
0)
Solve[{Eqn1,Eqn2,Eqn3,Eqn4,Eqn5},{r1,r2,b1,b2,p}]
===
Subject: Re: Finding unknown parameters using Mathematica
Hi Shafiq.
there is something wrong in your code. What should:=E2=88=91 mean? I
think Sum. Note that it is better if you change your code to InputForm
before posting it. This is done by selecting and then menu: Cell/Convert
To/InputForm.
Short Form of drivative == explicite Form of derivative
(D[L,r1]==...), what shlould give True. Note that the form: D[L,r1] is
only used in differential equations. In your case you would write:
eqns={ explicite form of derivative == 0, ...}
Solve[eqns,vars]
Daniel
===
Subject: Re: Limiting range of variables when defining funxtions of several
variables
The problem is what information is available when the test is
performed, in this case the test is performed when the pattern for y
is matched but x in the conditional expression is not the result of
pattern matching. In general if you need to test conditions that
depend on multiple arguments move the condition test to the right
hand side of the definition:
test[x_,y_]:=x+y/;y<=1-x
Ssezi
===
Subject: (Not trivial) Definite Integration of a rational function
I think it is time to open a new thread about integration, isn't it?
Anyway...Here we go..(Mathematica 5.2 is used).
f = HoldForm[Integrate[(2*x)/((x + 1)*(x^3 + 3*x^2 + 2*x + 1)), {x, 0,
Infinity}]]
We have
(*plot*)
machine precision*)
0.3712169752602472
Attempts to get an analytic result failed.
Timing[ReleaseHold[f]]
{145.969*Second, Integrate[(2*x)/((1 + x)*(1 + 2*x + 3*x^2 + x^3)),
{x, 0, Infinity}]}
{144.765*Second, Integrate[(2*x)/((1 + x)*(1 + 2*x + 3*x^2 + x^3)),
(
BTW, Mathematica version 4 returns Infinity for this definite
integral. I think this happens because a partial fraction
decomposition is performed and the integrand is written
Apart[(2*x)/((x + 1)*(x^3 + 3*x^2 + 2*x + 1))]
-(2/(1 + x)) + (2*(1 + 2*x + x^2))/(1 + 2*x + 3*x^2 + x^3)
and after it integrates definitely each term of last sum
(Integrate[#1, {x, 0, Infinity}] & ) /@ %
Integrate::idiv: Integral of 1/(1 + x) does not converge on
{0,Infinity}.
Integrate::idiv: Integral of (1 + x)^2/(1 + 2*x + 3*x^2 + x^3) does
not
converge on {0,Infinity}.
Integrate[-(2/(1 + x)), {x, 0, Infinity}] + Integrate[(2*(1 + 2*x +
x^2))/(1 + 2*x + 3*x^2 + x^3), {x, 0, Infinity}]
One additional curious thing is it doesn't show the relevant warning
message for convergence...
)
Let see the antiderivative
Together[D[%, x]]
Plot[%%, {x, 0, 10}]; (*plot*)
2*(-Log[1 + x] + RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[x - #1] +
2*Log[x - #1]*#1 + Log[x - #1]*#1^2)/
(2 + 6*#1 + 3*#1^2) & ])
(2*x)/((1 + x)*(1 + 2*x + 3*x^2 + x^3))
I wonder why Mathematica since it gets the indefinite integral, it
fails to evaluate
the definite integral. Any ideas?
Since both the integrand and the antiderivative posses no
singularities/discontinuities
in the integration range I proceed to get the definite integral by
application of the Newton-Leibniz
formula
This time we get finally the definite integral
Simplify[%[[2]]]
N[%]//Chop
{176.90699999999998*Second, -2*RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & ,
(Log[-#1] + 2*Log[-#1]*#1 + Log[-#1]*#1^2)/(2 + 6*#1 + 3*#1^2)
& ]}
-2*RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[-#1] + 2*Log[-#1]*#1 +
Log[-#1]*#1^2)/(2 + 6*#1 + 3*#1^2) & ]
0.37121697526024766
But timings like this look quite unreasonable for integrals like this
(at least for me!).
Do you have other ideas of getting this definite integral (with better
time performances
if possible)?
It is also looks quite strange that in two other CAS (the one can be
freely downloaded from Internet)
I got directly the definite integral in a couple of seconds.
Any explanation about previous time performance of Limit will be
greatly appreciate.
Trying to write down the result
-2*RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[-#1] + 2*Log[-#1]*#1 +
Log[-#1]*#1^2)/(2 + 6*#1 + 3*#1^2) & ]
in a expanded format (i e no RootSum object; of course I am aware that
it will be lost the compacteness
that RootSum offers) I tried
Tr[Simplify[(-2*((Log[-#1] + 2*Log[-#1]*#1 + Log[-#1]*#1^2)/(2 + 6*#1
+ 3*#1^2)) & ) /@
(x /. Solve[(1 + 2*#1 + 3*#1^2 + #1^3 & )[x] == 0, x])]]
Chop[N[%]]
-((2*(12 + 2^(1/3)*(27 - 3*Sqrt[69])^(2/3) + 6*3^(1/3)*(2/(9 -
Sqrt[69]))^(2/3))*
Log[1 + (2/(27 - 3*Sqrt[69]))^(1/3) + ((1/2)*(9 -
Sqrt[69]))^(1/3)/3^(2/3)])/
(3*(6 + 2^(1/3)*(27 - 3*Sqrt[69])^(2/3) + 6*3^(1/3)*(2/(9 -
Sqrt[69]))^(2/3)))) -
(2*(8*I*(9 - Sqrt[69])^(2/3) + 2^(1/3)*3^(1/6)*(2*2^(1/3)*3^(1/6)*(-
I + Sqrt[3]) +
(9 - Sqrt[69])^(1/3)*(-9 - 3*I*Sqrt[3] + I*Sqrt[23] +
Sqrt[69])))*
Log[1 + (-1 + I*Sqrt[3])/(2^(2/3)*(27 - 3*Sqrt[69])^(1/3)) - ((1 +
I*Sqrt[3])*((1/2)*(9 - Sqrt[69]))^(1/3))/(2*3^(2/3))])/
(3*(4*I*(9 - Sqrt[69])^(2/3) + 2^(1/3)*3^(1/6)*(2*2^(1/3)*3^(1/6)*(-
I + Sqrt[3]) +
(9 - Sqrt[69])^(1/3)*(-9 - 3*I*Sqrt[3] + I*Sqrt[23] +
Sqrt[69])))) -
(2*(8*(9 - Sqrt[69])^(2/3) + 2^(1/3)*3^(1/6)*(2*I*2^(1/3)*3^(1/6)*(I
+ Sqrt[3]) +
(9 - Sqrt[69])^(1/3)*(-9*I - 3*Sqrt[3] + Sqrt[23] +
I*Sqrt[69])))*
Log[1 + (-1 - I*Sqrt[3])/(2^(2/3)*(27 - 3*Sqrt[69])^(1/3)) + (I*(I
+ Sqrt[3])*((1/2)*(9 - Sqrt[69]))^(1/3))/(2*3^(2/3))])/
(3*(4*(9 - Sqrt[69])^(2/3) +
2^(1/3)*3^(1/6)*(2*I*2^(1/3)*3^(1/6)*(I + Sqrt[3]) +
(9 - Sqrt[69])^(1/3)*(-9*I - 3*Sqrt[3] + Sqrt[23] +
I*Sqrt[69]))))
0.3712169752602468
Am I missing something simpler than this setting?
Is it a built -in function that RootSum and this function acts like
the pair {RootReduce,ToRadicals} below?
Reduce[1 + 2*x + 3*x^2 + x^3 == 0, Reals]
ToRadicals[%]
RootReduce[%]
x == Root[1 + 2*#1 + 3*#1^2 + #1^3 & , 1]
x == -1 - (2/(3*(9 - Sqrt[69])))^(1/3) - ((1/2)*(9 - Sqrt[69]))^(1/3)/
3^(2/3)
x == Root[1 + 2*#1 + 3*#1^2 + #1^3 & , 1]
Dimitris
===
Subject: Re: Closing All Input Cells at Once- KB shortcuts
You can find the same thing in the help browser:
Perhaps this version is more printer-friendly.
===
Subject: Re: Closing All Input Cells at Once- KB shortcuts
Yes, it was the pdf I printed with the very nun-useful
results. I can use this document by eye but I can think of
several places it would be much more useful when printed.
I'm amazed that whoever made this document messed it up so