Can you pad the interior data points, by introducing duplicates which are only slightly off-value in the independent? That way it would technically be piecewise linear, but would almost completely behave like piecewise (near-)constant.

For example,

x:=[0.0,3.0,6.0,9.0]:
y:=[4.0,5.0,7.0,11.0]:

plots:-pointplot(zip(`[]`,x,y),style=line,view=0..11,scaling=constrained);

padx:=[x[1],seq([x[i]-1e-9,x[i]][],i=2..nops(x))]:
pady:=[seq([y[i],y[i]][],i=1..nops(x)-1),y[-1]]:

plots:-pointplot(zip(`[]`,padx,pady),style=line,view=0..11,scaling=constrained);

You might also make each new entry in pady be ever so slightly different that its (x-wise farther) neighbor, if say the piecewise linear interpolator doesn't like zero slope. This might not be necessary. Ie,

pady:=[seq([y[i],y[i]+1e-9][],i=1..nops(x)-1),y[-1]];

You could do,

plots:-listplot(result);

acer

That unknown _B1~ has the assumption that it should be either 0 or 1.

Substituting either of those two assumed values of _B1~ into your test gives something which simplifies to a trivial identity in both cases.

restart:

eq := sin(x^2)=1/2:
_EnvAllSolutions := true:

s := solve(eq);
                           1     (1/2)    1     (1/2)   
                      s := - (%1)     , - - (%1)        
                           6              6             
                    %1 := 24 Pi _B1~ + 72 Pi _Z1~ + 6 Pi

test := subs(x=s[1],eq);

                           /2                       1   \   1
                test := sin|- Pi _B1~ + 2 Pi _Z1~ + - Pi| = -
                           \3                       6   /   2

vars:=[ op( indets(test,name) minus {constants} ) ];

                            vars := [_B1~, _Z1~]

about(vars[1]);

  Originally _B1, renamed _B1~:
    is assumed to be: OrProp(0,1)


expand( subs(vars[1]=0, test) );

                                    1   1
                                    - = -
                                    2   2

expand( subs(vars[1]=1, test) );

                                    1   1
                                    - = -
                                    2   2

acer

If you don't upload the worksheet we can only guess, with so few details to go on.

There might be code in the Startup Code region.

There might be code in a collapsed Code Edit Region.

There might be action code behind an embedded component.

There might be code in a collapsed Execution Group within a Document Block.

And within or without any of those there might be calls to procedures or modules stored in separate .mla Library archive files.

acer

As Axel suggested, the imaginary component is zero.

Following that idea, one can simplify the expression (RHS-LHS of the equation), and justifiably consider its real component only. There appears to be a root where the expression's purely real value passes (continuously) through zero near dd = -2.98 or so.

restart;

TTot := 70:
TC := 17:
GM := .26:
QMax := 870:
V := 3600*GM*QMax*TTot:

Eq := V = Int(QMax*exp((-t+TC)/dd)*(1+(t-TC)/TC)^(TC/dd), t = 0 .. TTot):

# The imaginary part is zero
simplify(evalc(Im((rhs-lhs)(Eq))));

                               0

eq := value(Eq):
expr:=Re( simplify((convert((rhs-lhs)(eq),rational)),size) ):

H:=proc(x)
  if not type(x,numeric) then
    return 'procname'(args);
  end if;
  evalf(eval(expr,dd=x));
end proc:

Digits:=20:

sol := fsolve(Re(H(x)), x=-10..0);

                     -2.9847453623107653100

#plot(H(x),x=-9..-2);

evalf[1000](subs(dd=sol,Re(eq))):
evalf(%);

                        7                           7
          5.700240000 10  = 5.7002399999999999996 10 

The expression does not exceed about -5.7e7 for dd>0.

acer

Your expressions have multiple (at least two) different assumed names that print as p~.

indets(Q[6][2,3],name);

                           {p~, p~, x}

Those do not combine or cancel under `simplify` (which is usual).

I have not yet looked to see where they arise. It may be that you could use `additionally` rather than repeated calls to `assume`. Or you might be able to proceed with something crude like,

simplify(convert(Q[6][2,3],`global`));

[update]

You just have to move the line which calls `assume(p,real)` out of the P_chain procedure. It can be called just once at the top-level (since your P_chain already declares p as global), at before the loop. Here is an edited version,

June_18_modif.mw

acer

With Maple 18 one could try and use the new InertForm package for this.

The dot appearing in the imaginary term looks a little heavy to me. (It shows in the Std GUI, but not the TTY).

One nice thing is that the extra brackets don't appear, when done in the GUI.

For example,

  restart:

  c:=a+b*I:
  expr:=c^2:

  U := `%+`(evalc(Re(expr)),evalc(I*Im(expr))):

  with(InertForm):

  Display(U);

  Display(U, inert=false);

  Value(U);

acer

Are you using deprecated matrix rather than Matrix, and if so why?

Are you using Maple 18 and have you declared D as a top level local? If you are using name D at the top level (ie. this code is not appearing inside any procedure of your own) then perhaps it would simplify your issue if you switched that to say DD. 

How about something simpler such as this, if the Q[i] are all Matrices,

DD := add( evalf( (1/3)*Q[i]*(z(i+1)^3-z(i)^3) ), i=1..4 );

Note that even if you kept the loop, rather than use add, then you would not need evalm if you used Matrix instead of matrix.

If you still get problems with the diagonal elements, then perhaps you code upload the code so that we could investigate directly.

acer

I don't yet see a way to get an exact chatcterization of all the roots.

You may be able to obtain float approximations to roots. It may help to scale the expression (according to its magnitude on the requested range? Unclear to me.)

restart:

expr:=-9999990000000000000000*cos(166*beta)*sinh(166*beta)*cosh(88*beta)^2
 -9999990000000000000000*cos(88*beta)^2*sin(166*beta)*cosh(166*beta)
 +9999990000000000000000*sinh(166*beta)*cos(88*beta)^2*cos(166*beta)
 +9999990000000000000000*cosh(88*beta)^2*cosh(166*beta)*sin(166*beta)
 +10000010000000000000000*cos(166*beta)*sinh(166*beta)
 +10000010000000000000000*sin(166*beta)*cosh(166*beta)
 +9999990000000000000000*sinh(88*beta)*cos(166*beta)*cosh(166*beta)*cosh(88*beta)
 -9999990000000000000000*sinh(88*beta)*sin(166*beta)*sinh(166*beta)*cosh(88*beta)
 +9999990000000000000000*sin(88*beta)*cos(88*beta)*sinh(166*beta)*sin(166*beta)
 +9999990000000000000000*cos(88*beta)*cos(166*beta)*sin(88*beta)*cosh(166*beta)
 -9980010000000000000000*cosh(88*beta)^2*sinh(166*beta)*cos(88*beta)^2*cos(166*beta)
 -9980010000000000000000*cosh(88*beta)^2*cosh(166*beta)*sin(166*beta)
 *cos(88*beta)^2
 +9980010000000000000000*sinh(88*beta)*cos(88*beta)^2*sin(166*beta)
 *sinh(166*beta)*cosh(88*beta)
 -9980010000000000000000*cos(88*beta)*cosh(88*beta)^2*sin(88*beta)
 *sin(166*beta)*sinh(166*beta)
 +9980010000000000000000*sinh(88*beta)*cosh(88*beta)*cosh(166*beta)
 *cos(88*beta)^2*cos(166*beta)
 +9980010000000000000000*cosh(88*beta)^2*cosh(166*beta)*cos(88*beta)
 *sin(88*beta)*cos(166*beta)-9980010000000000000000*cos(88*beta)
 *sinh(88*beta)*cos(166*beta)*sin(88*beta)*sinh(166*beta)*cosh(88*beta)
 +9980010000000000000000*cos(88*beta)*cosh(88*beta)*sin(88*beta)
 *sin(166*beta)*cosh(166*beta)*sinh(88*beta):

#Digits:=500:
plot(expr,beta=0.0..0.2, view=0..1e-100);
#Digits:=10:

findroots:=proc(expr,a,b,{guard::posint:=5,maxtries::posint:=50})
local F,x,sols,i,res,start,t;
   x:=indets(expr,name) minus {constants};
   if nops(x)>1 then error "too many indeterminates"; end if;
   F:=subs(__F=unapply(expr,x[1]),__G=guard,proc(t)
      Digits:=Digits+__G;
      __F(t);
   end proc);
   sols,i,start:=table([]),0,a;
   to maxtries do
      i:=i+1;
      res:=RootFinding:-NextZero(F,start,
                                 'maxdistance'=b-start);
      if type(res,numeric) then
         sols[i]:=fnormal(res);
         if sols[i]=sols[i-1] then
            start:=sols[i]+1.0*10^(-Digits);
            i:=i-1;
         else
            start:=sols[i];
         end if;
      else
         break;
      end if;
   end do;
   op({entries(sols,'nolist')});
end proc:

findroots(expr, 0, 0.2); # missing some hinted at by the plot

  0.01641624893, 0.03039949207, 0.05572565081, 0.06691672633, 0.1327837617


findroots(expr*1e-90, 0, 0.2); # agreeing with the plot

 0.01641624893, 0.03039949207, 0.05572565081, 0.06691672633, 0.09465903392, 

   0.1039052938, 0.1327837505, 0.1327837510, 0.1417059480, 0.1700140953, 

   0.1804062466

acer

expr1 := p+p^(-1/(theta-1))*sum(q[i]^(theta/(theta-1)),i=(1..n));

                                         /  n                  \
                           /      1    \ |-----     /  theta  \|
                           |- ---------| | \        |---------||
                           \  theta - 1/ |  )       \theta - 1/|
             expr1 := p + p              | /    q[i]           |
                                         |-----                |
                                         \i = 1                /

expr2 := p^(-1/(theta-1))*(p^(theta/(theta-1))
         +sum(q[i]^(theta/(theta-1)),i=1..n));

                              /               /  n                  \\
                /      1    \ | /  theta  \   |-----     /  theta  \||
                |- ---------| | |---------|   | \        |---------|||
                \  theta - 1/ | \theta - 1/   |  )       \theta - 1/||
      expr2 := p              |p            + | /    q[i]           ||
                              |               |-----                ||
                              \               \i = 1                //

H:=(ee,a) -> a*simplify(1/a*ee):

H(expr1, p^(-1/(theta-1)));

                          /               /  n                  \\
            /      1    \ | /  theta  \   |-----     /  theta  \||
            |- ---------| | |---------|   | \        |---------|||
            \  theta - 1/ | \theta - 1/   |  )       \theta - 1/||
           p              |p            + | /    q[i]           ||
                          |               |-----                ||
                          \               \i = 1                //

Did you want to also programmaticaly "discover" that second argument to H above, the special term, perhaps as the coefficient of the Sum subexpression?

acer

Is this what you're after?

FirstList:=[A,B,C,D,E,F]:

combinat:-choose(FirstList,4);
   [[A, B, C, D], [A, B, C, E], [A, B, C, F], [A, B, D, E], 

     [A, B, D, F], [A, B, E, F], [A, C, D, E], [A, C, D, F], 

     [A, C, E, F], [A, D, E, F], [B, C, D, E], [B, C, D, F], 

     [B, C, E, F], [B, D, E, F], [C, D, E, F]]

nops(%);
                               15

acer

 Another possibility might be to use Maple's CodeGeneration[Matlab] command to produce code for a Matlab .m function from the Maple procedure `symb_sol` (then written out to a file), so that all the subsequent pieces are done right in the running Matlab session. 

We haven't been shown the equations which are solved ( via eliminate) so its not possible to tell yet whether all of `symb_sol` can be so handled.

acer

It's difficult to figure out the whole situation, as we do not seem to know what SystemEquations_NoDriver are.

But I don't under stand why you are forming all the symbolic input equations and holding BB_Driverlist_set(i) for all `i` in Maple memory at the same time.

You run (in Matlab), a first loop to store every symbolic equation set, in the BB_Driverlist_set(i). This is the step you've described as loading the control points into Maple.

Then you run another loop (in Matlab), to apply the Maple symb_sol function to every data set previously stored. Each time through the loop this produces a set of results as equations. And then you take rhs's and stuff the numeric results into the next part of a single Matrix. It's a bit dodgy about ordering of equations in a set, according to the name on the lhs, but it works.

Using a single Matrix is probably a good idea. How about a single loop which takes the input set (for a given time point), applies symb_sol, and then puts the result into the appropriate slot in the single result Matrix. And do not store or assign BB_Driverlist_set(i) for all i or for all SimNum.

And especially do not ever create any equation of the form name=expression, in a loop.

I also do not understand why you have to form the symbolic equations in sets. Perhaps I misunderstand the inputs. (Recall that we don't know what your SystemEquations_NoDriver is.)

As I understand it, you are calling from a running Matlab session into Maple to execute various steps of solving and manipulation. By "solving" you seem to mean that you generate an explicit general symbolic solution, and then evaluate this at each timestep's input data. (We'll have to take your word on this being more desriable than repeated numeric solving.)

It seems to me that you are doing something a bit like this,

pars:={ICE_T, BAT_P, EM2_T, EM2_W}:

eqs:={ICE_T=var1+cos(var2),BAT_P=sin(var1),EM2_T=var1^2,EM2_W=var3^2};

                 /                               2              2  
         eqs := { BAT_P = sin(var1), EM2_T = var1 , EM2_W = var3 , 
                 \                                                 

                                   \ 
           ICE_T = var1 + cos(var2) }
                                   / 

symb_sol:=unapply(eqs,[var1,var2,var3]):

ansMat := Matrix(nops(pars), 50000): #Matrix(m, numtimesteps):

# Loop over these next steps...

for i from 1 to 1 do

  ans := symb_sol(1.1,2.3,3.7);  print(ans);

  ansMat[..,i] := Vector([seq( rhs(ans[i]), i=1..nops(pars) )]):

end do:

  {BAT_P = 0.8912073601, EM2_T = 1.21, EM2_W = 13.69, ICE_T = 0.4337239787}

ansMat[..,1];

                               [0.8912073601]
                               [            ]
                               [        1.21]
                               [            ]
                               [       13.69]
                               [            ]
                               [0.4337239787]

How about making it a bit more like this,

pars:={ICE_T, BAT_P, EM2_T, EM2_W}:

eqs:={ICE_T=var1+cos(var2),BAT_P=sin(var1),EM2_T=var1^2,EM2_W=var3^2}:

exprs:=eval([BAT_P, EM2_T, EM2_W, ICE_T], eqs);

                     [               2      2                  ]
            exprs := [sin(var1), var1 , var3 , var1 + cos(var2)]

symb_sol:=unapply(exprs,[var1,var2,var3]):

ansMat := Matrix(nops(pars), 50000): #Matrix(m, numtimesteps):

# Loop over these next steps...

for i from 1 to 1 do

  ans := symb_sol(1.1,2.3,3.7);  print(ans);

  ansMat[..,i] := Vector(ans,datatype=float[8]):

end do:

                  [0.8912073601, 1.21, 13.69, 0.4337239787]

ansMat[..,1];

                             [0.891207360100000]
                             [                 ]
                             [ 1.21000000000000]
                             [                 ]
                             [ 13.6900000000000]
                             [                 ]
                             [0.433723978700000]

Those are just ideas for improvement. Better still might be to have no named equations, anywhere. Pass float[8] Vectors or Matrices of timestep data into Maple. Pump these directly into a procedure. (You may need help in altering your simple `unapply` into another kind of procedure...) And put the results directly into the slots of a float[8] answer Matrix. No equations, no calls to `rhs`, no sets, no lists.

How about passing a single float[8] Matrix of size numtimesteps-by-numinputs in one call from Matlab to Maple. Then make another call from Matlab to Maple which runs a Maple procedure (which you write) that evaluates the "solver" on each row of input data, and then places the results into a row of a results float[8] Matrix. Then pass the whole float[8] answer Matrix back to Maple. Ie, about one or three calls from Matlab to Maple, not 50000+.

acer

Note:

restart:             

'10.2^20';

                                              21
                               0.1485947396 10

Even special evaluation rules (to prevent premature evaluation of arguments in a procedure call, say) will not prevent automatic simplification. Having said that,

restart:

p := proc() 10.2^20; end proc;

                         p := proc() 10.2^20 end proc

evalf[50](p());
                                                             21
               0.14859473959783543420355740092833203224576 10


restart:       

evalf[50](`^`(10.2,20));

                                                             21
               0.14859473959783543420355740092833203224576 10

acer

Using the parameters and compile option of dsolve/numeric allows the same two last plots to be produced, with the overall computation time reduced from about 75 sec to 38 sec.

I also changed the 2D Math of the DE system to 1D Maple Notation. This seems to make the GUI use less resources (and respond better while editing). Your mileage may vary.

Benzene_mod.mw

acer