acer

32303 Reputation

29 Badges

19 years, 309 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

You can supply an option to the call to the plot command, for specifying the relative location of the first axis.

And you can change the size option below, to get dimensions more to your liking.

restart

k := 0.1e-1:

Pinit := 100:

V := 25:

rho := 1:

eq := diff(P(t), t) = -k*P(t)+rho*Cp*(int(V*sin(alpha), t));

diff(P(t), t) = -0.1e-1*P(t)+12.5*sin(5)*t

solution := dsolve({eq, P(0) = Pinit}, P(t));

P(t) = 1250*sin(5)*t-125000*sin(5)+exp(-(1/100)*t)*(100+125000*sin(5))

plot(rhs(solution), t = 0 .. 600, axis[1] = [location = low], size = [500, 400])

Download time-dependent_aerodynamic_pressure_ac.mw

See the Help page for topic  plot,options  for more details on 2D plotting choices.

The following reuable procedure P produces an expression (which is intended as pretty-printing in Maple's GUI similarly as does the the unevaluated ContinuedFraction function call).

Are you able to utilize the following expression assigned to res (either by copying its pretty-printed 2D output as special, rich-text, etc, or as generated MathML?)

restart

f := NumberTheory:-ContinuedFraction((3+sqrt(7))*(1/2)); lprint(%)

_m139948790579680

NumberTheory:-ContinuedFraction(3/2+1/2*7^(1/2))

P := proc (ee, N::posint) local f, H, i, T; H := Term(NumberTheory:-ContinuedFraction(ee), 0 .. N-1); T := `..`; for i from N by -1 to 2 do if (H[i])::list then T := H[i, 1]/(H[i, 2]+T) else T := 1/(H[i]+T) end if end do; H[1]+T end proc

res := P((3 + sqrt(7))/2, 10);

2+1/(1+1/(4+1/(1+1/(1+1/(1+1/(4+1/(1+1/(1+1/(1+`..`)))))))))

lprint(MathML:-ExportPresentation(res));

"<math xmlns='http://www.w3.org/1998/Math/MathML'><mrow><mn>2</mn><mo>+</mo><mf\
rac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><mfrac><mn>1</mn><mrow><mn>4</mn><mo>+<\
/mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><mfrac><mn>1</mn><mrow><mn>1</mn\
><mo>+</mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><mfrac><mn>1</mn><mrow><m\
n>4</mn><mo>+</mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><mfrac><mn>1</mn><\
mrow><mn>1</mn><mo>+</mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><mi>..</mi>\
</mrow></mfrac></mrow></mfrac></mrow></mfrac></mrow></mfrac></mrow></mfrac></mr\
ow></mfrac></mrow></mfrac></mrow></mfrac></mrow></mfrac></mrow></math>"

#res := P(cos(x), 12);

#ContinuedFraction(cos(x), x);

#MathML:-ExportPresentation(res);

 

Download ContFrac_MathML.mw

Your worksheet has some coding mistakes.

1) The definition of prof1 has some function calls where multiplcation was likely intended. Eg,
   (...)(ps - s)
is is a function call, and not a product of two brackets terms. For multiplication you'd need either an explicit multiplication symbol, or (in 2D Input) a space, between the adjacent closing&opening brackets )( .  This occurs a few times.

2) The attempted multiple assignments like,
   [p,q] := solve(...)
is invalid. This occurs a few times.

It's often useful to test your procedure with some actual float arguments, before trying to plot it. That would reveal to you that WhyNot3 wasn't working at all.

I've corrected those aspects, in this attachment.  HmaxProc_ac.mw

Here's another approach.

If you'd prefer the result as just a list of two numbers then you could leave off the [a,b]=~ bit at the start.

L := [a,b]=~applyrule(a::numeric*f^b::numeric=[a,-b],X1);

          L := [a = 1.234, b = 6.789]

A couple of other examples,

[a,b]=~applyrule(a::numeric*f^b::numeric=[a,-b],-1/(f^6.789));
                      [a = -1, b = 6.789]

[a,b]=~applyrule(a::numeric*f^b::numeric=[a,-b],f^6.789);
                      [a = 1, b = -6.789]

Here's your original example, round-tripped,

X := a/f^b;

a/f^b

X1 := 1.234/(f^6.789);

1.234/f^6.789

L := [a,b]=~applyrule(a::numeric*f^b::numeric=[a,-b],X1);

[a = 1.234, b = 6.789]

eval(X, L);

1.234/f^6.789

Download rt.mw

ps. The order of terms in the result is forced by us, here, since the order of terms in the returned list is specified in the code.

pps. You could also tweak the specified types, eg. only looking for negative b, etc. You could also wrap it up in a reusable procedure. You could also use indets (and a structured type) to yank out multiple such subexpressions from a larger expression. Hard to say more because we don't yet know exactly how you intend to utilize this.

@Samuel Hesselberg Any new document will have its own Startup code region (which is empty, at first). The Startup code is not automatically inherited from any other already opened sheet.

Your original .mw file is a document with its own, separate Startup code.

So if you copy&paste the main body of that older document then you'll also have to provide access to the code in the older document's Startup region as well. The embedded components you're pasting in rely upon and call that code.

You can add it to the Startup code region of your newly opened document, or you can add in somewhere else in your new document's body (any document block, or even a code-edit-region).

First, find the Startup code in your original. Then figure out how you want your new document to be able to also use it. The easiest way is also to copy the older document's Startup code to the new document's Startup code region.

I am not sure about a good approach for Maple 2015.

In Maple 2018 and later you could use plots:-textplot with its rotation option.

An example is given below.

Note that I did not have the lines extend down to y=0, since your image didn't either. I also didn't make it handle the lines if f(x) becomes negative. I figure that if you want such a more general example or different behaviour you would provide it -- and all your requirements -- explicitly and clearly.

restart;

with(plots):

f := x -> exp(cos(x)+sin(x));

proc (x) options operator, arrow; exp(cos(x)+sin(x)) end proc

(a,b) := 0,3:

n := 12:

h := (b-a)/n:

L := [seq([x,evalf(f(x))],x=a..b,h)]:

plots:-display(
  plot(f,a..b),
  pointplot(L[2..-2],symbolsize=10,symbol=solidcircle),
  seq(textplot([L[i][],
                typeset(" ",eval('Typesetting:-Typeset'('f'(evalf[3](L[i][1]))))
                            =evalf[5](L[i][2])),
                rotation=Pi/2],align=above,font=[DEFAULT,10]),i=2..nops(L)-1),
  seq(plottools:-line([L[i][1],min(L[..,2])],L[i],thickness=0),i=1..nops(L)-1),
  tickmarks=[L[..,1],[$ceil(min(L[..,2]))..floor(max(L[..,2]))]],
  font=[DEFAULT,9], size=[700,450]);

Download rotated_textplot_ex2.mw

Are you asking about codegen:-optimize? You haven't usefully provided an example of what exactly you're trying to accomplish.

If it's replacing generated names in results from codegen:-optimize then it's not clear to me whether you're dealing with Maple expressions of procedures. Also, while codegen:-optimize does have a documented mechanism to utilize a given list of preferred names, it's somewhat awkward IMO.

Perhaps something like the following might serve:

restart;

 

rep := (L::list(`=`),s::symbol)
  -> subs((nm->nm=parse(s||(String(nm)[2..])))~(lhs~(L)),
          L)[]:

 

repproc := (P::procedure,s::symbol)
  -> subs((nm->nm=parse(s||(String(nm)[2..])))~([op(2,eval(P))]),
          eval(P)):


We use rep to replace the temporary names in the optimized computation sequence.

f := x*ln(x)+2*x^2*ln(x)+3*x^3*ln(x);

x*ln(x)+2*x^2*ln(x)+3*x^3*ln(x)

temp := codegen:-optimize(f);

t1 = ln(x), t3 = x^2, t9 = 3*t1*t3*x+2*t1*t3+t1*x

rep([temp], b);

b1 = ln(x), b3 = x^2, b9 = 3*b1*b3*x+2*b1*b3+b1*x

 

We use repproc to replace the locals in the optimized procedure.

 

p := proc(x) x*ln(x)+2*x^2*ln(x)+3*x^3*ln(x); end proc;

proc (x) x*ln(x)+2*x^2*ln(x)+3*x^3*ln(x) end proc

temp2 := codegen:-optimize(p);

proc (x) local t1, t3; t1 := ln(x); t3 := x^2; 3*t1*t3*x+2*t3*t1+t1*x end proc

repproc(temp2, b);

proc (x) local b1, b3; b1 := ln(x); b3 := x^2; 3*b1*b3*x+2*b3*b1+b1*x end proc

Download cg_rep.mw

You have p(lambda__1 + p) which is a function call (to p), and not a multiplication. Perhaps you intended that as p*(lambda__1 + p) instead.

restart

psi = Pi^2*n^2+p; q := diff(theta(y), y, y) = psi^2*theta(y)

bc := theta(0) = 0, theta(1) = `&lambda;__1`*(cos(n*Pi)-1)/(p*(`&lambda;__1`+p)*n)

q := dsolve({bc, q})

theta(y) = lambda__1*(cos(n*Pi)-1)*exp(psi*y)/(n*p*(exp(psi)*p+exp(psi)*lambda__1-p*exp(-psi)-lambda__1*exp(-psi)))-lambda__1*(cos(n*Pi)-1)*exp(-psi*y)/(n*p*(exp(psi)*p+exp(psi)*lambda__1-p*exp(-psi)-lambda__1*exp(-psi)))

simplify(convert(q, trigh))

theta(y) = lambda__1*(cos(n*Pi)-1)*sinh(psi*y)/(p*(lambda__1+p)*n*sinh(psi))

Download nima_ac.mw

To prevent  -2*x-3*x from automatically simplifying to -5*x, etc, you can use an inert form.

Here, the inert form only needs to retain inert `%+`.

You can also force the visual display to have + in the usual black, instead of the gray that denotes the inert operator.

restart;

with(InertForm):

p := x^2 - a*x - b*x +a*b ;

a*b-a*x-b*x+x^2

L := [[2,3],[-3,7],[9,10]];

[[2, 3], [-3, 7], [9, 10]]

inexpr := eval(MakeInert(sort(p,x)),[`%^`=`^`,`%*`=`*`]);

`%+`(x^2, -a*x, -b*x, a*b)

temp := eval~(inexpr,map[2](Equate,[a,b],L))

[`%+`(x^2, -2*x, -3*x, 6), `%+`(x^2, 3*x, -7*x, -21), `%+`(x^2, -9*x, -10*x, 90)]

Display~(temp, inert=false);

0, "%1 is not a command in the %2 package", _Hold, Typesetting

value(temp);

[x^2-5*x+6, x^2-4*x-21, x^2-19*x+90]

Download inert_ex2.mw

If you add
   option remember;
to FirmModelFC and FirmModelPP then it won't waste as much time in duplicate computations, generating both the 3rd and 4th items of their returned expression sequences.

Also, you can guard against premature evaluation by having diffr1 and diffr2 return unevaluated when passed mere symbolic names (not yet numbers) as its arguments. That allows you to use the expression form calling sequence of plots:-inequal, as opposed to operator form calling sequence.

I removed the commented stuff, just to make it more legible to me.

I changed the range for delta to 0..1/2, since I myself didn't see too much of great interest for delta>1/2. Change that back if wanted, of course.

restart

with(plots)

``

c := 1; cr := 0.3e-1*c; u := 1; sExp := 0.6e-1*c; s := .65*c; v := 3*c

``

FirmModelPP := proc (alpha, delta) local p0, xi0, q0, Firmpf0, G0, Recpf0, Unsold0, Environ0; option remember; xi0 := 1; p0 := min(s+sqrt((v-s)*(c-s)), delta*v+sExp); q0 := u*(v-p0)/(v-s); f(N) := 1/u; F(N) := N/u; G0 := int(F(N), N = 0 .. q0); Firmpf0 := (p0-c)*q0-(p0-s)*G0; Recpf0 := (sExp-cr)*xi0*q0; Environ0 := q0+G0; Unsold0 := G0; return p0, q0, Firmpf0, Recpf0, Unsold0, Environ0 end proc

``

FirmModelFC := proc (alpha, beta, delta) local p00, xi00, q00, Firmpf00, G00, Recpf00, Unsold00, Environ00, pr00; option remember; xi00 := 1; p00 := s+sqrt((v-s)*(c-s)); if p00 < delta*v+sExp then q00 := u*(v-p00)/(v-s); f(N) := 1/u; F(N) := N/u; G00 := int(F(N), N = 0 .. q00); Firmpf00 := (p00-c)*q00-(p00-s)*G00; Recpf00 := `&xi;00*q00*`(sExp-cr); Unsold00 := G00; Environ00 := q00+Unsold00 else q00 := alpha*u*(v-p00)/(v-s); f(N) := 1/u; F(N) := N/u; G00 := int(F(N), N = 0 .. q00/alpha); pr00 := p00-delta*v; Firmpf00 := (p00-c)*q00-alpha*(p00-s)*G00; Recpf00 := (beta*(pr00-sExp)+sExp-cr)*xi00*q00-(1/2)*(pr00-sExp)*beta^2*xi00^2*q00^2/(u*(1-alpha)); Unsold00 := G00; Environ00 := q00+Unsold00 end if; return p00, q00, Firmpf00, Recpf00, Unsold00, Environ00 end proc

``

diffr1 := proc (alpha, delta) if not [alpha, delta]::(list(numeric)) then return ('procname')(args) end if; FirmModelPP(alpha, delta)[3]-FirmModelFC(alpha, .2, delta)[3] end proc

diffr2 := proc (alpha, delta) if not [alpha, delta]::(list(numeric)) then return ('procname')(args) end if; FirmModelPP(alpha, delta)[4]-FirmModelFC(alpha, .2, delta)[4] end proc

 

P1 := implicitplot(diffr1, 0 .. 1, 0 .. 1/2, color = black, thickness = 2)

NULL

P2 := implicitplot(diffr2, 0 .. 1, 0 .. 1/2, color = black, thickness = 2)

P3 := inequal({diffr1(alpha, delta) < 0, diffr2(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1/2, color = red)

NULL

P4 := inequal({diffr2(alpha, delta) > 0, diffr1(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1/2, color = yellow)

NULL

P5 := inequal({diffr1(alpha, delta) > 0, diffr2(alpha, delta) < 0}, alpha = 0 .. 1, delta = 0 .. 1/2, color = green)

NULL

P6 := inequal({diffr1(alpha, delta) > 0, diffr2(alpha, delta) > 0}, alpha = 0 .. 1, delta = 0 .. 1/2, color = "LightBlue")

NULL

display(P1, P2, P3, P4, P5, P6, scaling = constrained)

NULL

Download Code_JN_major_ac.mw

I didn't bother rewriting FirmModelPP and FirmModelAC altogether (eg. making evalhf'able counterparts, resolving the integration steps, etc). I suspect it could be faster still...

You appear to have used Maple 2021, and in that version the above all executed in about 9 seconds on my machine.

Do you mean something like the following?

I also added some sorting of the terms in the generated results, which I think makes the effect clearer.

note: Your given example might be handled by dealing merely with the second index j of the y[i,j] names. And the highest index accepted happens to match your example N values. But I deliberately put in the test of the first index, as well as sorting by the second index, so that it'd also handle a more general example (as you described it in words).

restart;

expr := r*y[0, 1]+y[0, 0]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/2)*r*(r-1)*y[-1, 2]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/720)*r*(r-1)*(r-2)*(r-3)*(r-4)*(1+r)*y[-3, 8]+(1/360)*r*(r-1)^2*(r-2)*(r-3)*(1+r)*y[-3, 7]+(1/720)*r*(r-1)*(r-2)*(r-3)*(1+r)*(2+r)*y[-3, 6];

r*y[0, 1]+y[0, 0]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/2)*r*(r-1)*y[-1, 2]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/720)*r*(r-1)*(r-2)*(r-3)*(r-4)*(1+r)*y[-3, 8]+(1/360)*r*(r-1)^2*(r-2)*(r-3)*(1+r)*y[-3, 7]+(1/720)*r*(r-1)*(r-2)*(r-3)*(1+r)*(2+r)*y[-3, 6]

 

T := proc(a,b)
       local anm:=indets(a,specindex(integer,y))[1],
             bnm:=indets(b,specindex(integer,y))[1];
       if op(1,anm)>op(1,bnm) then true;
       elif op(2,anm)<op(2,bnm) then true;
       else false; end if;
     end proc:

L := sort([op(expr)],T):

U := sort([indets(expr,specindex(integer,y))[]],T):

P := proc(n::nonnegint) local i;
       sort(add(L[i],i=1..n+1),order=plex(U[]));
     end proc:

 

P(0);

y[0, 0]

P(1);

y[0, 0]+r*y[0, 1]

P(2);

y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2]

P(3);

y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]

P(4);

y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]

P(5);

y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]

P(6);

y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]+(1/720)*r*(r-1)*(r-2)*(r-3)*(1+r)*(2+r)*y[-3, 6]

P(7);

y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]+(1/720)*r*(r-1)*(r-2)*(r-3)*(1+r)*(2+r)*y[-3, 6]+(1/360)*r*(r-1)^2*(r-2)*(r-3)*(1+r)*y[-3, 7]

P(8);

y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]+(1/720)*r*(r-1)*(r-2)*(r-3)*(1+r)*(2+r)*y[-3, 6]+(1/360)*r*(r-1)^2*(r-2)*(r-3)*(1+r)*y[-3, 7]+(1/720)*r*(r-1)*(r-2)*(r-3)*(r-4)*(1+r)*y[-3, 8]

 

Download sorty_thing.mw


Your given example could also be handled by the following simpler approach. But I think that the sorted effect is nicer. And if this is run in a separate session from the above code then the addends are not necessarily sorted so nicely (IMO), in each result.

expr := r*y[0, 1]+y[0, 0]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/2)*r*(r-1)*y[-1, 2]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/720)*r*(r-1)*(r-2)*(r-3)*(r-4)*(1+r)*y[-3, 8]+(1/360)*r*(r-1)^2*(r-2)*(r-3)*(1+r)*y[-3, 7]+(1/720)*r*(r-1)*(r-2)*(r-3)*(1+r)*(2+r)*y[-3, 6];

r*y[0, 1]+y[0, 0]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/2)*r*(r-1)*y[-1, 2]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/720)*r*(r-1)*(r-2)*(r-3)*(r-4)*(1+r)*y[-3, 8]+(1/360)*r*(r-1)^2*(r-2)*(r-3)*(1+r)*y[-3, 7]+(1/720)*r*(r-1)*(r-2)*(r-3)*(1+r)*(2+r)*y[-3, 6]

q := (n::nonnegint)->select(u->op(2,indets(u,indexed)[1])<=n,expr):

q(2);

r*y[0, 1]+y[0, 0]+(1/2)*r*(r-1)*y[-1, 2]

q(4);

r*y[0, 1]+y[0, 0]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/2)*r*(r-1)*y[-1, 2]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]

q(7);

r*y[0, 1]+y[0, 0]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/2)*r*(r-1)*y[-1, 2]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/360)*r*(r-1)^2*(r-2)*(r-3)*(1+r)*y[-3, 7]+(1/720)*r*(r-1)*(r-2)*(r-3)*(1+r)*(2+r)*y[-3, 6]

Download sorty_thing2.mw

You could get it with calls to GAMMA instead of WhittakerM.

f := exp(-beta*x)*x^(5-alpha):

simplify(convert(int(f,x),GAMMA));

((beta*x)^alpha*x^(-alpha)*(alpha-5)*(alpha-1)*(alpha-2)*(alpha-3)*(alpha-4)*GAMMA(-alpha+1, beta*x)-exp(-beta*x)*beta*(alpha-5)*(alpha-2)*(alpha-3)*(alpha-4)*x^(-alpha+1)+exp(-beta*x)*beta^2*(alpha-5)*(alpha-3)*(alpha-4)*x^(-alpha+2)-exp(-beta*x)*beta^3*(alpha-4)*(alpha-5)*x^(-alpha+3)+exp(-beta*x)*beta^4*(alpha-5)*x^(-alpha+4)-x^(5-alpha)*exp(-beta*x)*beta^5+x^(-alpha)*(beta*x)^alpha*GAMMA(-alpha+6))/beta^6

Download Wh_ex.mw

You could further collect on, say, exp, GAMMA, beta, etc. But mostly that just rearranges those terms. You haven't said what exactly you mean by less "complicated".

You can use the so-called 2-argument calling sequence of the eval command, to evaluate the expression at specific values for names upon which it depends.

If you are already have the lists of values then you can do it all programmatically, instead of typing out all those individual equations. That'd let you do it more easily if you have further collections of values to use.

expr := (`&varepsilon;__1`*exp(`&zeta;__1`*phi)+`&varepsilon;__2`*exp(zeta__2*phi))/(`&varepsilon;__3`*exp(`&zeta;__3`*phi)+`&varepsilon;__4`*exp(`&zeta;__4`*phi))

eval(expr, [`&varepsilon;__1` = 1, `&varepsilon;__2` = 1, `&varepsilon;__3` = 1, `&varepsilon;__4` = 0, `&zeta;__1` = 0, `&zeta;__2` = 1, `&zeta;__3` = -2, `&zeta;__4` = 0])

(1+exp(phi))/exp(-2*phi)

vars1 := [`&varepsilon;__1`, `&varepsilon;__2`, `&varepsilon;__3`, `&varepsilon;__4`]; vars2 := [`&zeta;__1`, `&zeta;__2`, `&zeta;__3`, `&zeta;__4`]

L1, L2 := [1, 1, 1, 0], [0, 1, -2, 0]

eval(expr, [Equate(vars1, L1)[], Equate(vars2, L2)[]])

(1+exp(phi))/exp(-2*phi)

Download eval_ex_1.mw

note. You could also construct a reusable procedure to do the substitutions, by using unapply on the expression (w.r.t. all those names). That doesn't really seem better to me, especially if you're only going to use it once.

expr := (`&varepsilon;__1`*exp(`&zeta;__1`*phi)+`&varepsilon;__2`*exp(zeta__2*phi))/(`&varepsilon;__3`*exp(`&zeta;__3`*phi)+`&varepsilon;__4`*exp(`&zeta;__4`*phi))

F := unapply(expr, `&varepsilon;__1`, `&varepsilon;__2`, `&varepsilon;__3`, `&varepsilon;__4`, `&zeta;__1`, `&zeta;__2`, `&zeta;__3`, `&zeta;__4`)

F(1, 1, 1, 0, 0, 1, -2, 0)

(1+exp(phi))/exp(-2*phi)

Download eval_ex_2.mw

Hopefully I haven't misunderstood all your goals.

I did the first attempt (with a finer resolution of points) in two ways: one like your original, and one more efficient that stops (for each X value) when it finds the first W value that succeeds.

At the end I also threw in some stuff with going straight for a curve (finding the first W value of equality, for the given range). And also a bit with some extended ranges. Might be useful; I'm not sure. I don't know what the end-point -1.7 might mean to you, as a restriction...

restart;

Digits:=15;

15

Lambda[0]:=2/(1*(5+3*w))

2/(5+3*w)

h:=1*a^(-3/2*(1+w))

a^(-3/2-(3/2)*w)

L:=simplify(int(1/h,a))

2*a^(5/2+(3/2)*w)/(5+3*w)

F:=1/(sqrt(2*Pi)*tau)*exp(-(lambda-Lambda[0])^2/(2*tau^2))

(1/2)*2^(1/2)*exp(-(1/2)*(lambda-2/(5+3*w))^2/tau^2)/(Pi^(1/2)*tau)

F2:=simplify(subs(lambda=L,F)*diff(h,a)*1/a)

-(3/4)*a^(-7/2-(3/2)*w)*2^(1/2)*(1+w)*exp(-2*(a^(5/2+(3/2)*w)-1)^2/((5+3*w)^2*tau^2))/(tau*Pi^(1/2))

F3:=isolate(Lambda[0]-tau=subs(a=A,L),tau)

tau = -2*A^(5/2+(3/2)*w)/(5+3*w)+2/(5+3*w)

F4:=Int(simplify(combine(subs(F3,F2))),a=0.01..10,method=_d01ajc)
     assuming A>0.5, A<1, a>0.01, a<10, w>-1.7, w<-1;

Int((3/8)*a^(-7/2-(3/2)*w)*2^(1/2)*(1+w)*exp(-(1/2)*(a^(5/2+(3/2)*w)-1)^2/(A^(5/2+(3/2)*w)-1)^2)*(5+3*w)/(Pi^(1/2)*(A^(5/2+(3/2)*w)-1)), a = 0.1e-1 .. 10, method = _d01ajc)

valid:= proc(n,b,w,A)
if n<b then #print("ok, n,b",n,b,w,A);
return true:
else
return false
end if:
end proc:

X:=[seq(A,A=0.5..0.99,0.01)]:

nops(X)

50

W:=[seq(w,w=-1.7..-1,0.01)]:

nops(W);

71

forget(evalf); forget(`evalf/int`):
st:=time[real]():
S1:=[seq(ListTools:-Search(true,[seq(valid(evalf(F4),1/(8*rhs(F3)^2),w,A),
                                                       w=W)]),
                           A=X)]:
S1:=eval(S1,0=NULL);
time[real]()-st;

[68, 68, 68, 67, 66, 65, 63, 57, 46, 41, 38, 35, 32, 29, 26, 23, 20, 17, 13, 10, 6, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

2.873

forget(evalf); forget(`evalf/int`):
st:=time[real]():
T := 'T':
for ii from 1 to nops(X) do
  for jj from 1 to nops(W) do
    #print(ii,jj,W[jj],X[ii]);
    if evalf(eval(F4<1/(8*rhs(F3)^2),[w=W[jj],A=X[ii]])) then
      T[ii] := jj;
      #print([X[ii],W[jj]]);
      # value found for this A, bail out of the inner loop
      jj:=nops(W); next;
    end if;
  end do;
end do;
S1alt:=convert(T,list);
time[real]()-st;

[68, 68, 68, 67, 66, 65, 63, 57, 46, 41, 38, 35, 32, 29, 26, 23, 20, 17, 13, 10, 6, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

.887

Y1:=[seq(W[S1[i]],i=1..nops(S1))]:

P1:=plots:-pointplot([X,Y1],view=[0.5..1,-2..0],style=pointline,
              symbol=solidcircle,symbolsize=10,color=red,size=[500,300])

Y1alt:=[seq(W[S1alt[i]],i=1..nops(S1alt))]:

P1:=plots:-pointplot([X,Y1alt],view=[0.5..1,-2..0],style=pointline,
              symbol=solidcircle,symbolsize=10,color=red,size=[500,300])

Digits:=10:
plots:-display(
  plot(-1.7,0.5..1,linestyle=dot,color=gray),
plot(Xval->RootFinding:-NextZero(unapply('Re@evalf'(eval(-F4+1/(8*rhs(F3)^2),A=Xval)),w),-2.7,
                                 maxdistance=2.7,guardDigits=1, signchange=false, abstol=1e-5),
     0.5 .. 0.8, adaptive=true, view=[0.5..1,-2..0], size=[500,300])
);

unapply('evalf'(eval(-F4+1/(8*rhs(F3)^2),A=0.95)),w):
plot([Re@%,Im@%],-5..-1, size=[500,300]);

 

Download NumericScale_ac.mw

By "Gauss-Jordan form" do you mean reduced row echelon form, as produced by Gauss-Jordan elimination?

If so, then consider computing the r.r.e.f. of the Matrix augmented by the identity-matrix. The row operations that take the Matrix to its r.r.e.f. would then also be done to the identity matrix.

Below, if B the r.r.e.f. of M turns out to be the identity-matrix, then A (the result of the same row operations, but done to the identity matrix) is the inverse of M.

See also the subsection "Finding the inverse of a matrix" of this page on Gaussian elimination.

restart;
randomize():
with(LinearAlgebra):

n := 3:
M := RandomMatrix(n,generator=-3..3);

Matrix(3, 3, {(1, 1) = -2, (1, 2) = -1, (1, 3) = 1, (2, 1) = -1, (2, 2) = -3, (2, 3) = 1, (3, 1) = -3, (3, 2) = -3, (3, 3) = -3})

(B,A):=(a->(a[..,..n],
            a[..,n+1..]))(ReducedRowEchelonForm(<M|IdentityMatrix(n)>));

B, A := Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1}), Matrix(3, 3, {(1, 1) = -1/2, (1, 2) = 1/4, (1, 3) = -1/12, (2, 1) = 1/4, (2, 2) = -3/8, (2, 3) = -1/24, (3, 1) = 1/4, (3, 2) = 1/8, (3, 3) = -5/24})

A.M - B;

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0})

M^(-1);

Matrix(3, 3, {(1, 1) = -1/2, (1, 2) = 1/4, (1, 3) = -1/12, (2, 1) = 1/4, (2, 2) = -3/8, (2, 3) = -1/24, (3, 1) = 1/4, (3, 2) = 1/8, (3, 3) = -5/24})

M := RandomMatrix(n,generator=-3..3);

Matrix(3, 3, {(1, 1) = -3, (1, 2) = -1, (1, 3) = 1, (2, 1) = 0, (2, 2) = -2, (2, 3) = -1, (3, 1) = -1, (3, 2) = 1, (3, 3) = 1})

(B,A):=(a->(a[..,..n],
            a[..,n+1..]))(ReducedRowEchelonForm(<M|IdentityMatrix(n)>));

B, A := Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = -1/2, (2, 1) = 0, (2, 2) = 1, (2, 3) = 1/2, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0}), Matrix(3, 3, {(1, 1) = 0, (1, 2) = -1/2, (1, 3) = -1, (2, 1) = 0, (2, 2) = -1/2, (2, 3) = 0, (3, 1) = 1, (3, 2) = -2, (3, 3) = -3})

A.M - B;

Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0})

M^(-1);

Error, (in rtable/Power) singular matrix

Download some_la.mw

First 38 39 40 41 42 43 44 Last Page 40 of 336