tomleslie

Try...

Answered: tomleslie 13876

December 18 2018

0 4

restart;
expr:= (L__b+lambda)^2*(I*sqrt(3)*lambda-2*L__b+lambda)^7*(I*sqrt(3)*lambda+2*L__b-lambda)^7
/
(16384*(L__b^2-lambda*L__b+lambda^2)^5*lambda^3);
simplify(expand(expr));

Resultant simplification is pretty impressive!

#############################
# How I'd (probably?) do it.
#
# Initialize and get version
#
  restart;
  with(plots):
  kernelopts(version);
#
  imp_fun := x^3+y^3 = 9*x*y;;
  s:= evalf
      ( solve
        ( eval
          ( imp_fun,
            x=2
          )
        )
      );
  m1:= evalf
       ( eval
         ( implicitdiff
           ( imp_fun, y, x
           ),
           [ x=2, y=s[1] ]
         )
       );
#
# Plot the curve and the solution points for x=2
#
  cols:=[ red, green, blue, black]:
  display
  ( [ implicitplot
      ( imp_fun,
        x=-5..5,
        y=-5..5,
        color=blue,
        gridrefine=4
      ),
     pointplot
     ( [ seq
         ( [2,s[j]],
           j =1..numelems([s])
         )
       ],
       color=[ red, green, blue],
       style=point,
       symbol=solidcircle,
       symbolsize=20
     )
   ]
);

>

Download impDiff2.mw

Observation...

Answered: tomleslie 13876

December 17 2018

0 0

Seems to work just fine in Maple 18?????

See the attached

>	############################# # Initialize and get version # restart; kernelopts(version); # # OP's code # imp_fun := -4x + 10(x^2)*(y^(-2)) + y^2 =11; c := 2; s:= evalf( solve( subs( x = c, imp_fun))); m1 := evalf( subs ( { x = c, y = s[1] }, implicitdiff( imp_fun, y,x)));

-4*x+10*x^2/y^2+y^2 = 11

(1)

>

#############################
# How I'd (probably?) do it.
#
# Initialize and get version
#
  restart;
  with(plots):
  kernelopts(version);
#
  imp_fun := -4*x + 10*(x^2)*(y^(-2)) + y^2 =11;
  s:= evalf
      ( solve
        ( eval
          ( imp_fun,
            x=2
          )
        )
      );
  m1:= evalf
       ( eval
         ( implicitdiff
           ( imp_fun, y, x
           ),
           [ x=2, y=s[1] ]
         )
       );
#
# Plot the curve and the solution points for x=2
#
  display
  ( [ implicitplot
      ( imp_fun,
        x=-5..5,
        y=-5..5,
        color=blue,
        gridrefine=4
      ),
      plot
      ( [ seq
          ( [2,j],
            j in s
          )
        ],
        style=point,
        color=red,
        symbol=solidcircle,
        symbolsize=20
      )
    ]
  );

-4*x+10*x^2/y^2+y^2 = 11

>

Download impDiff.mw

More of a couple of observations, really...

Answered: tomleslie 13876

December 16 2018

1 1

Code executes "correctly" (ie with no errors) in Maple 18, Maple 2015, Maple 2016, Maple2017 and Maple 2018 - see the attached. Which Maple version are you using
In all of the above versions, the determinant of the system to be solved is very small (in some loop iterations being returned as identically zero), which means that the system has no unique solution, and a parameterized solution is returned. See the attached.

>

$restart; with(LinearAlgebra); with(Student[MultivariateCalculus]); f := unapply(x1*exp(t*x2-y), x1, x2); tt := Matrix(5, 1, [1, 2, 4, 5, 8]); yy := Matrix(5, 1, [3.2939, 4.2699, 7.1749, 9.3008, 20.259]); ff := f(2.50, .25); for ii to 5 do t := tt[ii][1]; y := yy[ii][1]; Jf := Gradient(f(x1, x2), [x1, x2] = [2.50, .25]); JF := Transpose(Jf[1]); JR := Transpose(Transpose(ff).JF); HF := Transpose(JF).JF; det := Determinant(HF); P := LinearSolve(HF, -JR); printf("%a \n", P) end do$

[Vector(2, {(1) = 0.4764869674e-1, (2) = .1191217418})]

Vector(2, {(1) = 0.5675995750168781e-2, (2) = 0.14189989369465868e-1})

Matrix(2, 2, {(1, 1) = 0.2270398301020486e-2, (1, 2) = 0.5675995750168781e-2, (2, 1) = 0.5675995750168781e-2, (2, 2) = 0.14189989369465868e-1})

Vector(2, {(1) = 0., (2) = -1.0})

Vector(2, [0.,-1.], datatype = float[8])

[Vector(2, {(1) = 0.2305436861e-1, (2) = .1152718430})]

Vector(2, {(1) = 0.13287597798991115e-2, (2) = 0.6643798896613762e-2})

Matrix(2, 2, {(1, 1) = 0.5315039120057534e-3, (1, 2) = 0.26575195588760484e-2, (2, 1) = 0.26575195588760484e-2, (2, 2) = 0.1328759778861665e-1})

$Vector(2, {(1) = -2.499999999783121-4.999999997831213*_t0[1], (2) = _t0[1]})$

Vector(2, [(-2.49999999978312)-4.99999999783121*_t0[1],_t0[1]])

[Vector(2, {(1) = 0.2081014008e-2, (2) = 0.2081014008e-1})]

Vector(2, {(1) = 0.10826548253730561e-4, (2) = 0.10826548253730561e-3})

Matrix(2, 2, {(1, 1) = 0.4330619301492224e-5, (1, 2) = 0.43306193014922243e-4, (2, 1) = 0.43306193014922243e-4, (2, 2) = 0.43306193014922245e-3})

Vector(2, {(1) = 0., (2) = -.25})

Vector(2, [0.,-.25], datatype = float[8])

[Vector(2, {(1) = 0.3188467430e-3, (2) = 0.3985584288e-2})]

Vector(2, {(1) = 0.2541581138e-6, (2) = 0.3176976422951935e-5})

Matrix(2, 2, {(1, 1) = 0.1016632455e-6, (1, 2) = 0.1270790569180774e-5, (2, 1) = 0.1270790569180774e-5, (2, 2) = 0.1588488211675247e-4})

$Vector(2, {(1) = -2.5000000000000004-12.500000001568154*_t2[1], (2) = _t2[1]})$

Vector(2, [(-2.5)-12.5000000015682*_t2[1],_t2[1]])

[Vector(2, {(1) = 0.1175484901e-7, (2) = 0.2350969802e-6})]

Vector(2, {(1) = 0.3454411881e-15, (2) = 0.6908823761e-14})

Matrix(2, 2, {(1, 1) = 0.1381764752e-15, (1, 2) = 0.2763529505e-14, (2, 1) = 0.2763529505e-14, (2, 2) = 0.5527059010e-13})

Vector[column](%id = 18446744074600137662)

Vector(2, [(-2.49999999957464)-20.*_t3[1],_t3[1]])

>

Download GNmethod.mw

You should always...

Answered: tomleslie 13876

December 14 2018

2 7

post code using the big green up-arrow in the Mapleprimes tool bar, because retyping your code from a picture is lengthy and error-prone.

When I do retype your code, it executes with no problems, so either

there is something else in your worksheet which i causing a problem
I have "inadvertently" fixed something because I typed what I expected to be there, rather than just copying your code character-by-character

If you can't figure out why the attached works and your original doesn't, then post your original using the big green u-arrow in the toolbar

>

tol:=1e-06;
n:=10;
f:=x->x^2-3;
x[0]:=0.1;
h:=x->2*x;
for k from 1 to n do
    x[k]:=evalf(x[k-1]-f(x[k-1])/h(x[k-1]));
    if abs(x[k]-x[k-1])< tol
    then print("Number of iterations"=k);
         print("approximate solution"=x[k]);
         print( f(x[k]) );
         break;
    end if;
    x[k-1]:=x[k];
end do:
plot( [f(x), f(x[k-1])+(x-x[k-1])*h(x[k-1])], x=-1..3, colour=[red,black]);

>

Download newtMeth.mw

Your fundamental problem...

Answered: tomleslie 13876

December 13 2018

0 1

is that on the first pass through your innermost loop 'tj' is an unassigned variable which you solve for.

On your second pass through the innermost loop 'tj' is no longer unassigned (because you solved for it on the first pass!). Now you can get around this by adding by adding tj:='tj' as the first statement in the loop, ensuring that on each loop iteration 'tj' is unassigned. If you do this, then all nested loops will execute with no problems, but the value of the variable 'som' will remain at zero, becuase the 'if' conditions which allow it to be updated are never satisfied - evidenced by adding a slight amendment to the print statement you had at this point ( which never executes)

I'm guessing that this is probably not the answer you want??? But see the attached anyway. NB this takes about 90secs to run on my machine - which I guess is fair enough to evaluate the innermost loop 36000 times. Something tells me there must be a better way: any more than three nested loops and I get a nosebleed!

>

restart;
som:=0:
for b1 from 10 to 10 by 1 do
  for b2 from 1 to 2 by 1 do
    for alpha from 0.5 to 0.5 by 0.1 do
      for beta from 0.33 to 0.5 by 0.1 do
        for c from 1 to 1 by 1 do
          for f from 1 to 10 by 1 do
            for g from 1 to 10 by 0.1 do
              for lambdaj from 0.2 to 0.4 by 0.1 do
                for gammaj from 0.2 to 0.4 by 0.1 do
                    p:='p';
                    tj:='tj'; # Added this
                    aiSQ:=(alpha*b1)/(alpha*b2+beta*b2+c);
                    ajSQ:=(beta*b1)/(alpha*b2+beta*b2+c);
                    UiSQ:=(1/2)*alpha*b1^2*(alpha^2*b2+2*alpha*beta*b2+c*alpha+beta^2*b2+2*beta*c)/(alpha*b2+beta*b2+c)^2;
                    UjSQ:=(1/2)*beta*b1^2*(alpha^2*b2+2*alpha*beta*b2+2*c*alpha+beta^2*b2+beta*c)/(alpha*b2+beta*b2+c)^2;
                    USQ:=(1/2)*b1^2*(alpha+beta)*(alpha*b2+beta*b2+2*c)/(alpha*b2+beta*b2+c)^2;
                    UTSQ:=UiSQ+UjSQ+USQ;
                    ai:=(-alpha*b2*f*p+alpha*b1*c-b2*f*p+b1*c)/(c*(alpha*b2+b2*beta+b2+c));
                    aj:=(alpha*b2*f*p+b1*beta*c+b2*f*p+c*f*p)/(c*(alpha*b2+b2*beta+b2+c));
                    aineg:=-(alpha*b2*f*lambdaj*p+b2*f*lambdaj*p+alpha*b1*c+b1*c)/(c*(alpha*b2*lambdaj-2*alpha*b2-b2*beta+b2*lambdaj-2*b2-c));
                    ajneg:=(alpha*b2*f*lambdaj*p+alpha*b1*c*lambdaj+b2*f*lambdaj*p+c*f*lambdaj*p-alpha*b1*c-b1*beta*c+b1*c*lambdaj-b1*c)/(c*(alpha*b2*lambdaj-2*alpha*b2-b2*beta+b2*lambdaj-2*b2-c));
                    uj:=beta*(b1*(aineg+ajneg)-(b2/2)*(aineg+ajneg)^2)-(c/2)*ajneg^2-p*f*(ajneg-aj);
                    uL:=(alpha+1)*(b1*(aineg+ajneg)-(b2/2)*(aineg+ajneg)^2)-(c/2)*aineg^2+p*f*(ajneg-aj);
                    eqtj:=gammaj*(uL-USQ)-((1-gammaj)/(1-lambdaj))*(uj-UjSQ)-tj;
                    tj:=solve(eqtj,tj);
                    dai:=diff(ai,p);
                    daj:=diff(aj,p);
                    daineg:=diff(aineg,p);
                    dajneg:=diff(ajneg,p);
                    dtj:=diff(tj,p);
                    Ujp:=beta*(b1*(ai+aj)-(b2/2)*(ai+aj)^2)-(c/2)*aj^2-p*f*(ajneg-aj)+(1-lambdaj)*tj;
                    ULp:=(alpha+1)*(b1*(ai+aj)-(b2/2)*(ai+aj)^2)-(c/2)*ai^2+p*f*(ajneg-aj)-tj-g*((p^2)/2);
                    eqp:=diff(ULp,p);
                    eqpp:=diff(eqp,p);
                    p:=solve(eqp,p);
                    CSQ:=b1-b2*(aiSQ+ajSQ);
                    Cabat:=b1-b2*(ai+aj);
                    Cneg:=b1-b2*(aineg+ajneg);
                    if   (ai>aineg)
                    then f*(aineg-ai)=0
                    else if   (aj>ajneg)
                         then f*(ajneg-aj)=0
                         else if (CSQ>0 and Cabat>0 and Cneg>0 and eqpp<0 and p>0 and p<1 and beta<alpha and aiSQ>0 and ajSQ>0 and ai>0 and aj>0 and aineg>0 and ajneg>0 and tj>0)
                              then #
                                   # Modifies the following print statement - which
                                   # never executes
                                   #
                                    printf("Made it to here %a %a %a %a %a %a %a %a %a %a\n",
                                            b1,b2,alpha,beta,c,f,g,lambdaj,gammaj,p
                                          );
                                   som:=som+1;
                              fi;
                         fi;
                    fi;
                od;
              od;
            od;
          od;
        od;
      od;
    od;
  od;
od;
som;

(1)

>

Download multiloop.mw

Anything wrong with...

Answered: tomleslie 13876

December 13 2018

0 2

int(cos(2*x)/(1+2*sin(3*x)^2), x):
eval(%, x=Pi)-eval(%, x=0)

which returns 0, more or less instantaneously

As written...

Answered: tomleslie 13876

December 13 2018

0 0

sin(x)*cos(x)

is a 2D plot? Did you possibly mean sin(x)*cos(y) ?

Any plot can be exported as EPS. The easiest way to do this, is

Generate the plot
Select and display the appropriate context menu by right clicking on the plot
Select "Encapsulated Postscript" from the "Export" sub-menu: follow the dialog

Suggest...

Answered: tomleslie 13876

December 12 2018

2 1

you examine the attached

>

NULL

>

Matrix(%id = 18446744074370527830)

(1)

>

Matrix(%id = 18446744074370524574)

(2)

>

Matrix(%id = 18446744074370521198)

(3)

>

Matrix(%id = 18446744074370507710)

(4)

>

expr := A.`<,>`(1, 2, 0) = `<,>`(0, -1, 2); eqns := [`θ__2` >= 0, `θ__2` <= Pi, `θ__3` >= 0, `θ__3` <= Pi, seq(lhs(expr)[j] = rhs(expr)[j], j = 1 .. 3)]; ans := solve(eqns); (eval(A, ans)).`<,>`(0, -1, 2)

eqns := [0 <= `θ__2`, `θ__2` <= Pi, 0 <= `θ__3`, `θ__3` <= Pi, -cos(`θ__2`)+2*sin(`θ__2`)*sin(`θ__3`) = 0, -2*cos(`θ__3`) = -1, sin(`θ__2`)+2*cos(`θ__2`)*sin(`θ__3`) = 2]

$ans := {`θ__2` = (1/6)*Pi, `θ__3` = (1/3)*Pi}$

Vector[column](%id = 18446744074370493366)

(5)

>

Download transVec.mw

Just for amusement...

Answered: tomleslie 13876

December 12 2018

0 1

it is possible to do pretty much what you want (although I have no idea what information is conveyed)

>

  restart;
  with(GraphTheory):
#
# Enter "physical" position for graph vertices
#
  vp := [ [2.5, 21], [6, 13.5], [8, 10], [11, 24.5], [14.3, 19.4],
          [16.8, 26], [22, 21.5], [22, 17], [22.2, 12.5], [26.8, 23],
          [28, 20.5], [30, 25.5], [32, 21], [29.5, 16]
        ]:
#
# Initialise the Graph with the correct number
# of vertices
#
  G1:= Graph
       ( [ seq( j,
                j=1..numelems(vp)
              )
         ]
       ):
#
# Draw Graph with given "physical" positions
#
  SetVertexPositions(G1, vp);
  DrawGraph(G1);
#
# Define start node, then span the graph. As a
# by-product, return the total physical distance
# traversed
#
  startNode:=5:
  doSpan:= proc( sp::integer, phyD::list, dist:=0.0 )
                 uses Student:-Calculus1:
                 global G1:
                 local currPos,
                       dd,
                       pD:= phyD,
                       d:= (p, q)-> evalf( Distance(p, q) );
                 if   numelems(phyD)=2
                 then AddEdge( G1,
                               { pD[1][1], pD[2][1] }
                             );
                      return dist+d( pD[1][2], pD[2][2] ):
                 else currPos, pD:= selectremove(i->i[1]=sp, pD);
                      dd:= sort
                           ( [ seq
                               ( [ currPos[][1],
                                   j[1],
                                   d(currPos[][2], j[2])
                                 ],
                                 j in pD
                               )
                             ],
                             (x,y)-> evalb(x[3]<y[3])
                           )[1];
                      AddEdge( G1,
                               { dd[1], dd[2] }
                             );
                      doSpan( dd[2], pD, dist+dd[3] );
                 fi;
           end proc:

  doSpan( startNode,
          [ seq
            ( [j, vp[j]],
              j in Vertices(G1)
            )
          ]
        );
#
# Draw the connected graph
#
  SetVertexPositions(G1, vp);
  DrawGraph(G1);

>

Download physGraph.mw

Integrand contain two variables...

Answered: tomleslie 13876

December 12 2018

0 3

You can't numerically evaluate int(f(x,y)) with respect to y with x unknown.

For your specific case it is actually possible to get a (rather complicated) expression for int(f(x,y), x=0..infinity) - ie a symbolic representation containing the variable 'x', and then evaluate this for any(?) value of x, as in the attached

>

restart:

>

f1 := (-(1.671130827*10^(-21)*I)*((0.6132469529e-1-.5253537767*I)*x^3+(.2894208344-2.479398027*I)*y*x^2+(0.1107017247e-1+(-.2377238011-0.2774956673e-1*I)*y^2)*x+(-.1167302837+.9999999999*I)*y^3+(0.2912028092e-1-.2494663781*I)*y)*exp(-1.732070784*10^(-15)*x*y)*((-0.6132469529e-1-.5253537767*I)*x^3+(.3115611798+2.669068985*I)*y*x^2+(-0.1107017247e-1+(-.3699321954+0.4318229002e-1*I)*y^2)*x+(.1455017102+1.246477826*I)*y^3+(0.2512356422e-1+.2152274756*I)*y)*exp(1.732070784*10^(-15)*x*y)*exp(-2.178669712*10^(-15)*y^2-3.534838336*10^(-16)*x^2+(1.000000000*10^(-9)*I)*x/y)/x)*y;

(1)

>	ans:=unapply(int(f1,y=0..infinity),x): ans(1); ans(10);

(2)

>

Download doInt.mw

If I understand correctly...

Answered: tomleslie 13876

December 12 2018

1 0

then the attached ought to do it

Download mormalize.mw

First problem...

Answered: tomleslie 13876

December 11 2018

0 1

Code you supplied contains the line

phin:=simplify(unapply(Transpose(psi).(hatA+A0),t)):

but 'psi' is undefined, hence so is Transpose(psi), and thus phin etc. This means that your final plot() command produces an error, because the function your are trying to plot contains the undefined variable 'psi'

On my machine it takes about 15seconds, before Maple "Errors". Is this excessive?

Your immediate problem...

Answered: tomleslie 13876

December 11 2018

0 0

is with the first remove() statement in the procedure. A remove() statement expects its first argument to be a logical test which returns the value 'true' or 'false' for each entry in the second argument.. If the logical test returns 'true' then the corresponding entry in the second argument is removed: So for example

remove( i-> i>3, [1,2,3,4,5,6]);

will return [1,2,3]

I have added a print statement and a couple of comments in the attached to illustrate this issue. There are other simple(?) syntax issues eg (non-exhaustive list)

There is a second remove() statement in your procedure wihich will fail for simlar reasons to the first
The Disatnce command takes three arguments - (graph, vertex, vertex)
Probably others, which I didn't spot during a quick read

>

Primmetje := proc (aantal, posities, begin)
                  local knopenover, knopeninmst, huidig, V, kaart, e, a;
                #
                # Inserted the following, just in case any
                # any of the commands from the GraphTheory
                # package, *want* to execute
                #
                  uses GraphTheory;
                  knopenover := [seq(i, i = 1 .. aantal)];
                  knopeninmst := {};
                  huidig := [0, 0];
                  if   begin <> {}
                  then V := [begin];
                       knopeninmst := knopeninmst union {V}
                  end if;
                #
                # Inserted print statement to illustrate
                # issue - first argument, ie 'V' is *not*
                # a logical test which returns 'true' or
                # 'false'
                #
                  printf("%a, %a\n", V, knopenover):
                  remove(V, knopenover); # <- This fails
                  kaart := Graph(aantal);
                  SetVertexPositions(kaart, posities);
                  while nops(knopeninmst) < aantal do
                        for e in knopeninmst do
                            for a in knopenover do
                                if   huidig = [0, 0]
                                              or
                                              Distance(posities[e], posities[a]) < Distance(posities[huidig[1]], posities[huidig[2]])# <-Number of Arguments!
                                then huidig:= [e, a];
                                              knopeninmst := knopeninmst union {a};
                                              remove(a, knopenover); # <- I'm betting this will fail
                                              AddEdge(kaart, huidig)
                                end if
                            end do
                        end do
                  end do
              end proc:
vp := [2.5, 21], [6, 13.5], [8, 10], [11, 24.5], [14.3, 19.4], [16.8, 26], [22, 21.5], [22, 17], [22.2, 12.5], [26.8, 23], [28, 20.5], [30, 25.5], [32, 21], [29.5, 16];
Primmetje(14, vp, 1);

[[6, 13.5]], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]

Error, (in Primmetje) invalid boolean expression: [[6, 13.5]]

>	remove( i-> i>3, [1,2,3,4,5,6]);

(1)

>

Download graphProb.mw

###################################################################################################

However I think(?) you have a couple of bigger "conceptual" problems waiting for you - particularly with your understanding of the GraphTheory:-SetVertexPositions() and the GraphTheory:-Distance() function.

The former is just so that you can control how a "Graph" is drawn. The latter measures "distance" between two vertices as the least number of 'edges' which have to be traversed in the 'arc' between the two vertices.It has nothing to do with the physical distance between these vertices in some drawn representation of the Graph.

So for example if you have a graph with vertex 'a' connected to vertex 'b', and vertex 'b' connected to vertex 'c', (and no other edges) then the "Distance" between 'a' and 'c' will be always be two, since you have to traverse two edges. It does not matter how you decide to "draw" this graph, position the vertices, whatever

From where you are now...

Answered: tomleslie 13876

December 10 2018

0 5

the easy way is just to stick everything in a loop and send the complete results from LSSolve() and SolveEquations() to separate Maple tables. Data can then be extracted from these tables to construct matricse of results, which can then be exported to Excel, using ExcelTools:-Export().

Two things to note

Worksheet for all 20 n-values takes about 5 minutes to run - so wait
Almost anything can be extracted from the table and reformatted/prettified/whatever whilst creting the final matrices. I just extracted the parameter name=value equations. But one could creatt the parameter names as column headers in the final Matrix with values as column entries - depends very much on what kind of format you want

>

(1)

>

Loading LinearAlgebra

>

(2)

>

(3)

>

(4)

>

(5)

>

$for n to 20 do e1 := p[s] = (n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3]))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n; e2 := p[b] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n; e3 := r[0] = lambda[0]*(20*p[0, p]*(1-r[0])*((7/2)*p[0, b]-(7/2)*p[0, b]*p[0, c]^5+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-30*p[0, c]^5)+20*r[0]*(7/2+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-(67/2)*p[0, c]^5)+p[0, b]*(11258.065*p[s]/p[b]+469.725)*(p[0, p]*(1-r[0])*((7/2)*p[0, b]-(7/2)*p[0, b]*p[0, c]^5+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-30*p[0, c]^5)+r[0]*(7/2+(15/2)*p[0, c]+(15/2)*p[0, c]^2+(15/2)*p[0, c]^3+(15/2)*p[0, c]^4-(67/2)*p[0, c]^5))+(469.725*(-p[0, p]*r[0]+p[0, p]+r[0]))*(-4*p[0, c]^5+p[0, c]^4+p[0, c]^3+p[0, c]^2+p[0, c])+11727.79-11727.79*p[0, c]^5*(-p[0, p]*r[0]+p[0, p]+r[0])+p[0, c]^5*(20*p[0, p]*(1-r[0])*(7*p[0, b]*(1/2)+30)+3018.625*r[0]+(20*p[0, p]*(1-r[0])*(7*p[0, b]*(1/2)+30)+670*r[0])*p[0, b]*(11258.065*p[s]+469.725*p[b])/(20*p[b])+2348.625*p[0, p]*(1-r[0]))); e4 := r[1] = lambda[1]*(20*p[1, p]*(1-r[1])*((15/2)*p[1, b]-(15/2)*p[1, b]*p[1, c]^5+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-62*p[1, c]^5)+20*r[1]*(15/2+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-(139/2)*p[1, c]^5)+p[1, b]*(11258.065*p[s]/p[b]+469.725)*(p[1, p]*(1-r[1])*((15/2)*p[1, b]-(15/2)*p[1, b]*p[1, c]^5+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-62*p[1, c]^5)+r[1]*(15/2+(31/2)*p[1, c]+(31/2)*p[1, c]^2+(31/2)*p[1, c]^3+(31/2)*p[1, c]^4-(139/2)*p[1, c]^5))+(469.725*(-p[1, p]*r[1]+p[1, p]+r[1]))*(-4*p[1, c]^5+p[1, c]^4+p[1, c]^3+p[1, c]^2+p[1, c])+11727.79-11727.79*p[1, c]^5*(-p[1, p]*r[1]+p[1, p]+r[1])+p[1, c]^5*(20*p[1, p]*(1-r[1])*(15*p[1, b]*(1/2)+62)+3738.625*r[1]+(20*p[1, p]*(1-r[1])*(15*p[1, b]*(1/2)+62)+1390*r[1])*p[1, b]*(11258.065*p[s]+469.725*p[b])/(20*p[b])+2348.625*p[1, p]*(1-r[1]))); e5 := r[2] = lambda[2]*(20*p[2, p]*(1-r[2])*((31/2)*p[2, b]-(31/2)*p[2, b]*p[2, c]^5+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-478*p[2, c]^5)+20*r[2]*(31/2+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-(987/2)*p[2, c]^5)+p[2, b]*(11258.065*p[s]/p[b]+469.725)*(p[2, p]*(1-r[2])*((31/2)*p[2, b]-(31/2)*p[2, b]*p[2, c]^5+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-478*p[2, c]^5)+r[2]*(31/2+(63/2)*p[2, c]+(127/2)*p[2, c]^2+(255/2)*p[2, c]^3+(511/2)*p[2, c]^4-(987/2)*p[2, c]^5))+(469.725*(-p[2, p]*r[2]+p[2, p]+r[2]))*(-4*p[2, c]^5+p[2, c]^4+p[2, c]^3+p[2, c]^2+p[2, c])+11727.79-11727.79*p[2, c]^5*(-p[2, p]*r[2]+p[2, p]+r[2])+p[2, c]^5*(20*p[2, p]*(1-r[2])*(31*p[2, b]*(1/2)+478)+12218.625*r[2]+(20*p[2, p]*(1-r[2])*(31*p[2, b]*(1/2)+478)+9870*r[2])*p[2, b]*(11258.065*p[s]+469.725*p[b])/(20*p[b])+2348.625*p[2, p]*(1-r[2]))); e6 := r[3] = lambda[3]*(20*p[3, p]*(1-r[3])*((31/2)*p[3, b]-(31/2)*p[3, b]*p[3, c]^5+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-478*p[3, c]^5)+20*r[3]*(31/2+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-(987/2)*p[3, c]^5)+p[3, b]*(11258.065*p[s]/p[b]+469.725)*(p[3, p]*(1-r[3])*((31/2)*p[3, b]-(31/2)*p[3, b]*p[3, c]^5+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-478*p[3, c]^5)+r[3]*(31/2+(63/2)*p[3, c]+(127/2)*p[3, c]^2+(255/2)*p[3, c]^3+(511/2)*p[3, c]^4-(987/2)*p[3, c]^5))+(469.725*(-p[3, p]*r[3]+p[3, p]+r[3]))*(-4*p[3, c]^5+p[3, c]^4+p[3, c]^3+p[3, c]^2+p[3, c])+11727.79-11727.79*p[3, c]^5*(-p[3, p]*r[3]+p[3, p]+r[3])+p[3, c]^5*(20*p[3, p]*(1-r[3])*(31*p[3, b]*(1/2)+478)+12218.625*r[3]+(20*p[3, p]*(1-r[3])*(31*p[3, b]*(1/2)+478)+9870*r[3])*p[3, b]*(11258.065*p[s]+469.725*p[b])/(20*p[b])+2348.625*p[3, p]*(1-r[3]))); e7 := t[0] = (-p[0, c]^5+1)/((1-p[0, c])*(r[0]+(1-r[0])/p[0, p]+(-p[0, c]^5+1)/(1-p[0, c])+(7*(p[0, b]*(1-r[0])+r[0]))/(2-2*p[0, b])+(15*p[0, c]^4+15*p[0, c]^3+15*p[0, c]^2+15*p[0, c])/(2-2*p[0, b]))); e8 := t[1] = (-p[1, c]^5+1)/((1-p[1, c])*(r[1]+(1-r[1])/p[1, p]+(-p[1, c]^5+1)/(1-p[1, c])+(15*(p[1, b]*(1-r[1])+r[1]))/(2-2*p[1, b])+(31*p[1, c]^4+31*p[1, c]^3+31*p[1, c]^2+31*p[1, c])/(2-2*p[1, b]))); e9 := t[2] = (-p[2, c]^5+1)/((1-p[2, c])*(r[2]+(1-r[2])/p[2, p]+(-p[2, c]^5+1)/(1-p[2, c])+(31*(p[2, b]*(1-r[2])+r[2]))/(2-2*p[2, b])+(511*p[2, c]^4+255*p[2, c]^3+127*p[2, c]^2+63*p[2, c])/(2-2*p[2, b]))); e10 := t[3] = (-p[3, c]^5+1)/((1-p[3, c])*(r[3]+(1-r[3])/p[3, p]+(-p[3, c]^5+1)/(1-p[3, c])+(31*(p[3, b]*(1-r[3])+r[3]))/(2-2*p[3, b])+(511*p[3, c]^4+255*p[3, c]^3+127*p[3, c]^2+63*p[3, c])/(2-2*p[3, b]))); e11 := p[0, p] = 1-exp(-lambda[0]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725)); e12 := p[1, p] = 1-exp(-lambda[1]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725)); e13 := p[2, p] = 1-exp(-lambda[2]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725)); e14 := p[3, p] = 1-exp(-lambda[3]*(-449.725*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+(11258.065*(n*t[0]/(1-t[0])+n*t[1]/(1-t[1])+n*t[2]/(1-t[2])+n*t[3]/(1-t[3])))*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n+469.725)); e15 := p[0, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[0]); e16 := p[1, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[1]); e17 := p[2, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[2]); e18 := p[3, c] = 1-((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[3]); e19 := p[0, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[0]))^1; e20 := p[1, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[1]))^1; e21 := p[2, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[2]))^2; e22 := p[3, b] = 1-(((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[3]))^6; e23 := p[0, s] = n*t[0]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[0]); e24 := p[1, s] = n*t[1]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[1]); e25 := p[2, s] = n*t[2]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[2]); e26 := p[3, s] = n*t[3]*((1-t[0])*(1-t[1])*(1-t[2])*(1-t[3]))^n/(1-t[3]); CSTR := `~`[`=`]({p[b], p[s], p[0, b], p[0, c], p[0, p], p[0, s], p[1, b], p[1, c], p[1, p], p[1, s], p[2, b], p[2, c], p[2, p], p[2, s], p[3, b], p[3, c], p[3, p], p[3, s], r[0], r[1], r[2], r[3], t[0], t[1], t[2], t[3]}, 0 .. 1); residuals := `~`[lhs-rhs]([e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12, e13, e14, e15, e16, e17, e18, e19, e20, e21, e22, e23, e24, e25, e26]); ansLS[n] := Optimization:-LSSolve(residuals, op(CSTR), iterationlimit = 10000); ansDS[n] := DirectSearch:-SolveEquations(residuals, CSTR, initialpoint = ansLS[n][2]) end do$