Advanced Search

acer

32333 Reputation

29 Badges

19 years, 318 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

Contact acer

MaplePrimes Activity


These are answers submitted by acer

ifactors...

acer 32333
September 16 2020
2 8

You should be able to use the results from the ifactors command, programmatically. (You can trivially massage its output into similar formats of your own invention.)

ifactors(725);

        [1, [[5, 2], [29, 1]]]

ifactors(1125);

         [1, [[3, 2], [5, 3]]]

ifactors(2048);

             [1, [[2, 11]]]

yes...

acer 32333
September 16 2020
1 4

If I understand your sentences rightly then yes, the inner procedure that is defined in the body of the outer procedure will be recreated each time the outer is invoked/called, and this is unnecessary overhead (computational effort and memory management) that could be avoided.

If instead the procedure A were simply defined at a higher level than procedure B, then you should be able to call A from within B, under the lexical scoping rules.

You could also implement both A and B within the same module, if you would prefer some compartmentalization.

You could also make the outer procedure be the local ModuleApply of a module B, where A is another local procedure within that module. That would make B into a so-called appliable module, where the name B could be called with arguments.

Would you like examples of all this? Have you read the Programming Guide (esp. the section on modules?)

ZeroMatrix...

acer 32333
September 16 2020
2 1

If you call the ZeroMatrix command without specifying its option

    compact=false

then (by default) the Matrix will be constructed with a special indexing function.

The effect of that is that the Matrix can only contain zeroes. (This is great for efficiency in some programming situations, since very large zero-matrices can take up little/constant memory.)

Much simpler is just to use the Matrix command instead of the ZeroMatrix command. The default fill-values of a call to Matrix are zero, and those don't have to be specified explicitly.

This is a guess, since you did not provide you code. It's almost always better to provide full code to reproduce. (The green up arrow in the Mapleprimes editor allows you to upload and attach documents. )

> 

restart;
> 

MZ1 := LinearAlgebra:-ZeroMatrix(3, 4);

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

> 

MZ1[2,3] := 17;

Error, attempt to assign non-zero to zero Matrix entry

> 

MZ2 := LinearAlgebra:-ZeroMatrix(3, 4, compact=false);

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

> 

MZ2[2,3] := 17;

17

> 

MZ2;

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

> 

MZ3 := Matrix(3, 4);

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

> 

MZ3[2,3] := 17;

17

> 

MZ3;

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

> 

 

If I set the "rtable" (ie. Matrix) associated with the above DataTable to be MZ1 then I'll get a popup error "Unable to change value" whenever I attempt to edit an entry and hit the Return/Enter key,

 

But if I associate MZ2 or MZ3 with the DataTable instead then I can edit the entries.

 

Creating MZ3 as above is much simpler than creating MZ2 as above.

 

Download zeromatrix_edit.mw

2D, 3D, etc...

acer 32333
September 15 2020
1 0

For the 2D case here are two short versions. (They could be shorter still, but I wanted to show a few of the intermediate assignments.)

First, one using operators.

with(plots):
F:=(x,y)->4*x^2+9*y^2;
G:=[D[1],D[2]](F)(2,1);
display(arrow([2,1],G),
        implicitplot(F(x,y)=F(2,1),x=-4..4,y=-2..2),
        scaling=constrained);

And one using expression form.

with(plots):
f:=4*x^2+9*y^2;
fval:=eval(f,[x=2,y=1]);
gval:=eval(zip(diff,f,[x,y]),[x=2,y=1]);
display(arrow([2,1],gval),
        implicitplot(f=fval,x=-4..4,y=-2..2),
        scaling=constrained);

Both of those can be adjusted using various options for the arrow dimensions, etc.

And 3D versions can also be made. The 2D implicit plot can be raised to a level height on the 3D surface, or a single contour line can be specified at the specific height. See attached.

gradient_arrow_operator.mw

gradient_arrow.mw

something...

acer 32333
September 15 2020
1 1

I think that your example is interesting, since the "size" of the desired target expression is smaller than that of the original expression, using any of these as norm:
   length,  `simplify/size/size`,  MmaTranslator:-Mma:-LeafCount
So it seems reasonable to ask why the "smaller" form is not found directly by some flavour of call to simplify.

My compressor oracle found these two ways:

> 

ex1:=(2*sqrt(N__s) + N__s + 1)/(2*N__s + 2);

(2*N__s^(1/2)+N__s+1)/(2*N__s+2)

> 

eval(collect(collect(ex1,[N__s],__P),[__P]),
     __P=(u->u));

1/2+N__s^(1/2)/(N__s+1)

> 

thaw(collect(subsindets(ex1,radical,freeze),
     [freeze~(indets(ex1,radical))[]],simplify));

1/2+N__s^(1/2)/(N__s+1)

Download simp_radical_ex.mw

You probably already realize this, but you should not use the symbolic option of the simplify command if you care about results branch cuts and always being correct. However, your current example and desired target are equivalent. I just wanted to mention it, since you made it sound as if simplify(...,symbolic) might be generally applicable/valid.

something...

acer 32333
September 15 2020
1 1

I took reciprocals (twice, naturally) in order to get op(2,X) and v into the exact same form. But the rest was not so bad.

It's always interesting when simplify (or the is command) requires extra help. [I will submit a Software Change Request, citing the examples.)

> 

restart

 

In Paul J. Nahin's book "The story of the square of " sqrt(-1) ", "Scipione del Ferro discovered the roots of the depressed cubic f (below, with p and q positive) by replacing x with u + v and by assuming that g (below) = 0.

del Ferro then found the values of u and v below

 

> 

f := x^3+p*x = q
> 

g := 3*u*v+p
> 

u := ((1/2)*q+((1/4)*q^2+(1/27)*p^3)^(1/2))^(1/3); v := -(-(1/2)*q+((1/4)*q^2+(1/27)*p^3)^(1/2))^(1/3)


Although Maple's roots of f (Sol below) do not resemble del Ferro's, its two ops satisfy g and evaluate to x = 2 when p and q are 6 and 20 respectively, as in the example Nahin gives in his book.

 

How is this so when Maple's roots seem so different from del Ferro's?

 

> 

Sol := solve([f, p > 0, q > 0], x)

piecewise(And(0 < p, 0 < q), [{x = (1/6)*(108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)-2*p/(108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)}], [])

> 

X := `assuming`([eval(x, Sol[1])], [p > 0, q > 0])

(1/6)*(108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)-2*p/(108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)

> 

`assuming`([simplify(combine(g))], [p > 0, q > 0])

0

> 

simplify(eval(X, [p = 6, q = 20]))

2

> 

`assuming`([simplify(combine(rationalize(X-u-v)))], [p > 0, q > 0])

0

> 

simplify(eval(X-u-v, [p = 6, q = 20]))

0

> 

`assuming`([simplify(combine(rationalize(op(1, X)-u)))], [p > 0, q > 0])

0

> 

`assuming`([simplify(combine(rationalize(op(2, X)-v)))], [p > 0, q > 0])

0

> 

simplify(rationalize(u), radical); op(1, X); evalb(% = `%%`)

(1/6)*(108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)

(1/6)*(108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)

true

> 

`assuming`([simplify(rationalize(v))], [p > 0, q > 0]); `assuming`([1/simplify(1/v, radical)], [p > 0, q > 0]); evalb(% = `%%`)

-(1/6)*(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)

-(1/6)*(-108*q+12*(12*p^3+81*q^2)^(1/2))^(1/3)

true

> 

`assuming`([is(combine(rationalize(op(1, X)-u)) = 0)], [p > 0, q > 0])

true

> 

`assuming`([is(combine(rationalize(op(2, X)-v)) = 0)], [p > 0, q > 0])

true

> 

NULL

Download CubicRootsacc.mw

one way...

acer 32333
September 15 2020
4 0

The line,
    combine(factor(expand(ee)))
will handle your expression, assigned below to ee.

If you have similar examples you might wish to restrict the combine by its exp option. If you have other examples where you don't want a blanket expand then you could target the original exp calls by using subsindets.

> 

restart;
> 

ee := exp(4*y) + 2*exp(2*y) + 1;

exp(4*y)+2*exp(2*y)+1

> 

factor(expand(ee))

((exp(y))^2+1)^2

> 

combine(factor(expand(ee)));

(exp(2*y)+1)^2

> 

combine(factor(expand(ee)),exp);

(exp(2*y)+1)^2

> 

 

Download exp_factor.mw

typeset...

acer 32333
September 14 2020
1 6

Is it because you misspelled "typeset"?

 

procedures and expressions...

acer 32333
September 14 2020
1 1

It's important to understand the basic kinds of structures and operations for differentiating.

Once you understand that you could move on to typeset "2D Math" representations. But IMO it's very important to understand what might be going on underneath.

note. You don't have to define a procedure such as f below, in order to obtain an expression such as
   y*sin(z)+2*cos(y)
You could also just type that expression out, assign that to some name, differentiate using the diff command, etc.

After you understand the following, you might look at these two older posts, here and here. (Note that I have not checked how those behave in modern versions.) There may also be support syntax for your input form from the Physics package. But you didn't attach your actual attempts.

> 

restart;
> 

# A procedure that takes two arguments (ie. has two
# parameters.

f := (y,z) -> y*sin(z)+2*cos(y);

proc (y, z) options operator, arrow; y*sin(z)+2*cos(y) end proc

> 

# An expression

f(y,z);

y*sin(z)+2*cos(y)

> 

# Differentiating the procedure with respect to its first
# parameter, by applying the differential operator `D`.

D[1](f);

proc (y, z) options operator, arrow; sin(z)-2*sin(y) end proc

> 

# Using D as previously, and then evaluating at the
# arguments (Pi/2,Pi).

D[1](f)(Pi/2,Pi);

-2

> 

# Differentiating the expression, with respect to y

diff(f(y,z),y);

sin(z)-2*sin(y)

> 

# Differentiating the expression, with respect to y,
# and then evaluating at the point y=Pi/2,z=Pi.

eval(diff(f(y,z),y), [y=Pi/2,z=Pi]);

-2

 

And we can repeat all the above, but this time differentiating with
respect to the second parameter of procedure f, or with respect to
the name `z` for the expression f(y,z).

> 

D[2](f);

proc (y, z) options operator, arrow; y*cos(z) end proc

> 

D[2](f)(Pi/2,Pi);

-(1/2)*Pi

> 

diff(f(y,z),z);

y*cos(z)

> 

eval(diff(f(y,z),z), [y=Pi/2,z=Pi]);

-(1/2)*Pi

 

In 2D Input mode, the following syntax will also create a procedure.

note(!) This syntax will not accomplish the same thing in 1D plaintext
input. In that mode you should use the syntax at the start of this worksheet.

 

> 

restart
> 

"f(y,z):=y*sin(z)+2*cos(y)"

proc (y, z) options operator, arrow, function_assign; y*sin(z)+2*cos(y) end proc

> 

f := (y,z) -> y*sin(z)+2*cos(y);

proc (y, z) options operator, arrow; y*sin(z)+2*cos(y) end proc

> 

 

Download procs_expressions_differentiation.mw

conformal...

acer 32333
September 14 2020
1 3

Firstly, you can compare your technique's results with those of the plots:-conformal command (using suitable parameters, of course).  See attached.

Next, we can cobble together a version of multiple calls to plots:-conformal, so that we attempt a shading by hue which is similar to his scheme. (One could also scale the intensity, the "V" of "HSV", but I have not done that.)

Also -- and I expect this may be key for you -- one can compare his scheme on log(z) versus yours for exp(z), and vice-versa, and his for arcsin(z) against your for sin(z). Do that explain the discrepancy to your satisfaction, as a matter of direction of the operation? (I do not know why he chose to display the inverse mapping...)

> 

restart;
> 

with(plots):
> 

xMax := 1: yMax := 1:
N := 10; step := abs(xMax)/N:
GL := proc(x, y) options operator, arrow; x+I*y end proc:
f := exp;
G := {}: for k from -N+1 to N+1 do
  G := `union`(G, {complexplot(f(GL(x, (k-1)*step)), x = -xMax .. xMax, color = brown)});
  G := `union`(G, {complexplot(f(GL((k-1)*step, y)), y = -yMax .. yMax, color = brown)})
end do:

10

exp

> 

P1:=display(G,scaling=constrained,size=[300,300]):
P1;

> 

P2:=plots:-conformal(f(z),z=-xMax-yMax*I..xMax+yMax*I,
                 grid=[2*N+1,2*N+1],color=brown,
                 scaling=constrained,size=[300,300]):
P2;

> 

# compare with   http://davidbau.com/conformal/#log(z)
#

(m,n):=16,16;
P3:=plots:-display(seq(seq(plots:-conformal(f(z),
               z=   (-m*1/16+(i-1)*(1/16)) + (-n*1/16+(j-1)*(1/16))*I
                 .. (-m*1/16+(i)*(1/16))   + (-n*1/16+(j)*(1/16))*I,
                 grid=[2,2],numxy=[2,2],
                 color=ColorTools:-Color("HSV",[(Pi/2+arctan(((i-1)-m)/m,
                                                             ((j-1)-n)/n))/(2*Pi),
                                                1,1]),
                 scaling=constrained,size=[300,300]),
    i=1..2*m+1),j=1..2*n+1)):
P3;

16, 16

> 

# alternate bookkeeping
# compare with  http://davidbau.com/conformal/#sin(z)
#
f := arcsin;
(a,b,c,d):=-1,-1,1,1;
(m,n):=(c-a)*16, (d-b)*16;
plots:-display(seq(seq(plots:-conformal(f(z),
               z=(a+(i-1)*(c-a)/(m-1))+(b+(j-1)*(d-b)/(n-1))*I
                 .. (a+(i)*(c-a)/(m-1))+(b+(j)*(d-b)/(n-1))*I,
                 grid=[2,2],numxy=[2,2],
                 color=ColorTools:-Color("HSV",[(Pi/2
                                                 +arctan((i-1-(m-1)/2)/(m-1),
                                                         (j-1-(n-1)/2)/(n-1)))/(2*Pi),
                                                1,1]),
                 scaling=constrained,size=[300,300]),
    i=1..m-1),j=1..n-1))

arcsin

-1, -1, 1, 1

32, 32

> 

# compare with  http://davidbau.com/conformal/#arcsin(z)
#
f := sin;
(a,b,c,d):=-1,-1,1,1;
(m,n):=(c-a)*16, (d-b)*16;
plots:-display(seq(seq(plots:-conformal(f(z),
               z=(a+(i-1)*(c-a)/(m-1))+(b+(j-1)*(d-b)/(n-1))*I
                 .. (a+(i)*(c-a)/(m-1))+(b+(j)*(d-b)/(n-1))*I,
                 grid=[2,2],numxy=[2,2],
                 color=ColorTools:-Color("HSV",[(Pi/2
                                                 +arctan((i-1-(m-1)/2)/(m-1),
                                                         (j-1-(n-1)/2)/(n-1)))/(2*Pi),
                                                1,1]),
                 scaling=constrained,size=[300,300]),
    i=1..m-1),j=1..n-1))

sin

-1, -1, 1, 1

32, 32

> 

# compare with  http://davidbau.com/conformal/#exp(z)
#
f := ln;
(a,b,c,d):=-1,-1,1,1;
(m,n):=(c-a)*16, (d-b)*16;
plots:-display(seq(seq(plots:-conformal(f(z),
               z=(a+(i-1)*(c-a)/(m-1))+(b+(j-1)*(d-b)/(n-1))*I
                 .. (a+(i)*(c-a)/(m-1))+(b+(j)*(d-b)/(n-1))*I,
                 grid=[2,2],numxy=[2,2],
                 color=ColorTools:-Color("HSV",[(Pi/2
                                                 +arctan((i-1-(m-1)/2)/(m-1),
                                                         (j-1-(n-1)/2)/(n-1)))/(2*Pi),
                                                1,1]),
                 scaling=constrained,size=[300,300]),
    i=1..m-1),j=1..n-1))

ln

-1, -1, 1, 1

32, 32

> 

# compare with  http://davidbau.com/conformal/#log(z)
#
f := exp;
(a,b,c,d):=-1,-1,1,1;
(m,n):=(c-a)*16, (d-b)*16;
plots:-display(seq(seq(plots:-conformal(f(z),
               z=(a+(i-1)*(c-a)/(m-1))+(b+(j-1)*(d-b)/(n-1))*I
                 .. (a+(i)*(c-a)/(m-1))+(b+(j)*(d-b)/(n-1))*I,
                 grid=[2,2],numxy=[2,2],
                 color=ColorTools:-Color("HSV",[(Pi/2
                                                 +arctan((i-1-(m-1)/2)/(m-1),
                                                         (j-1-(n-1)/2)/(n-1)))/(2*Pi),
                                                1,1]),
                 scaling=constrained,size=[300,300]),
    i=1..m-1),j=1..n-1))

exp

-1, -1, 1, 1

32, 32

> 

# compare with  http://davidbau.com/conformal/#arctan(z)
#

f := tan;
(a,b,c,d):=-1,-1,1,1;
(m,n):=(c-a)*16, (d-b)*16;
plots:-display(seq(seq(plots:-conformal(f(z),
               z=(a+(i-1)*(c-a)/(m-1))+(b+(j-1)*(d-b)/(n-1))*I
                 .. (a+(i)*(c-a)/(m-1))+(b+(j)*(d-b)/(n-1))*I,
                 grid=[2,2],numxy=[2,2],
                 color=ColorTools:-Color("HSV",[(Pi/2
                                                 +arctan((i-1-(m-1)/2)/(m-1),
                                                         (j-1-(n-1)/2)/(n-1)))/(2*Pi),
                                                1,1]),
                 scaling=constrained,size=[300,300]),
    i=1..m-1),j=1..n-1))

tan

-1, -1, 1, 1

32, 32

> 

 

Download conformalfun.mw

ps. If I were to attempt such a coloring scheme from scratch I'd try and do it more gracefully and efficiently, rather than clumping all these plots:-conformal calls with shared edges (and duplicated calculations). I did it this way partly because I was trying a crude mimicing of his "1/16" block, with your own technique also in my mind, and partly because I wanted to see the end result sooner and with less coding effort.

nn is posint...

acer 32333
September 13 2020
0 1

It seems as if you are trying to find the discrete minimum, for nn a positive integer in 1..5.

If I have interpreted the goal properly, then I suspect that the following might not surprise anyone (where each xx[ii] attaints its lower bound, for each posint nn).

Hence my use of n-n^2 in the first computation below (where even that min and seq call are avoidable if we realize it is decreasing).

Your original formulation is an abuse of the notation and syntax, IMHO. But below is (just) one way to deal with the programming.

> 

restart;
> 

min(seq(n-n^2, n=1..5));

-20

> 

W := proc(nn)
  local ii, rf;
  if not nn::posint then return 'procname'(args); end if;
  rf := add(xx[ii]^2-nn,ii=1..nn);
  return rf;
end proc:
> 

K := proc(nn::posint)
  Optimization:-Minimize(W(nn),seq(xx[ii]=1..10,ii=1..nn));
end proc:
> 

min(seq(K(i)[1], i=1..5));

-20.

> 

K(5);

[-20., [xx[1] = HFloat(1.0), xx[2] = HFloat(1.0), xx[3] = HFloat(1.0), xx[4] = HFloat(1.0), xx[5] = HFloat(1.0)]]

> 

 

Download discr_opt.mw

like this?...

acer 32333
September 13 2020
2 4

You could try something like this,

> 

restart;
> 

replace_all_C:=proc(expr::anything)
  subsindets(expr, ':-suffixed(_C, posint)',
           u->c[parse(convert(u,string)[3..])]);
end proc:
> 

sol:=dsolve(diff(x(t),t$3) = x(t));

x(t) = _C1*exp(t)+_C2*exp(-(1/2)*t)*sin((1/2)*3^(1/2)*t)+_C3*exp(-(1/2)*t)*cos((1/2)*3^(1/2)*t)

> 

tmp:=replace_all_C(sol)

x(t) = c[1]*exp(t)+c[2]*exp(-(1/2)*t)*sin((1/2)*3^(1/2)*t)+c[3]*exp(-(1/2)*t)*cos((1/2)*3^(1/2)*t)

> 

 

Download q_acc.mw

references...

acer 32333
September 13 2020
0 2

There are several interesting ways to plot the logistic map, and it's not completely clear from your wording what kind of plot you are after (though I might guess that you might be looking for a straightforward solution, as a beginner in Maple).

On this forum, see here, here, here, and here.

You can also see the Bifurcation command, in more recent Maple versions.

See also here, for the cobweb plot example near the bottom.

faq...

acer 32333
September 12 2020
0 1

This is a fequently asked question.

f:=x->x^2+x^3:

g:=unapply(diff(f(x),x),x);

   g := x -> 3*x^2+2*x

g:=D(f);
   g := x -> 3*x^2+2*x

fit...

acer 32333
September 11 2020
2 4

Let me know if this is the kind of thing that you wanted to accomplish.

I used the DirectSearch ver.2 package for the nonlinear fitting (since the NonlinearFit command does local optimization and without a tight set of parameter ranges it can converge to a local optimum that is a measurably worse fit). You can install the DirectSearch ver.2 package in the Maple 2020 GUI from the menubar's Maple Cloud login, or obtain it from the Maple Cloud website here, or for much older Maple versions install from the Application Center here.
 

> 

restart;
> 

with(LinearAlgebra):
> 

#data:=ExcelTools:-Import("D:/Maple/donnees_cycles_sin_lejeunes.xlsx"):
data:=ExcelTools:-Import(cat(kernelopts(homedir),"/mapleprimes",
                             "/donnees_cycles_sin_lejeunes.xlsx")):
> 

PK1Expe:=Column(data,1):
> 

fHz:=3: #fréquence
lambdaS:=0.5:
lambdaD:=0.3:
lambdaTrig:= piecewise(t < 1/fHz, 1+lambdaS*t*fHz, t >= 1/fHz, 1+lambdaS+lambdaD*sin((2*Pi*(t-1/fHz)*fHz))):

with(LinearAlgebra):
tensF:=Matrix([[lambdaTrig,0,0],[0,1/sqrt(lambdaTrig),0],[0,0,1/sqrt(lambdaTrig)]]):
tensdF:=diff(tensF,t):
tensFDev:=tensF-(1/3)*Trace(tensF)*Matrix(3,shape=identity):
tensB:=Multiply(tensF,Transpose(tensF)):
tensL:=Multiply(tensdF,MatrixInverse(tensF)):
tensD:=(1/2)*(tensL+Transpose(tensL)):
tensBe:=Matrix([[Be11(t),0,0],[0,1/sqrt(Be11(t)),0],[0,0,1/sqrt(Be11(t))]]):
tensBeDev:=tensBe-(1/3)*Trace(tensBe)*Matrix(3,shape=identity):

ode:=tensdBe[1,1]=-(2/3)*tensBe[1,1]*Trace(tensD)+Multiply(tensL[1,1],tensBe[1,1])+Multiply(tensBe[1,1],Transpose(tensL[1,1]))-(2/eta)*a*Multiply(tensBeDev[1,1],tensBe[1,1]): ivp:=[ode,Be11(0)=1]:
> 

tensdBe:=diff(tensBe,t):
> 

ode_sol:=dsolve(ivp,numeric,parameters=[a,eta]):
> 

Be11Num:=proc(aValue,etaValue,tValue)
   global __old_a, __old_eta;
   local res;
   if not [aValue,etaValue,tValue]::list(numeric) then
      return 'procname'(args);
   end if;
   if __old_a<>A_value and __old_eta<>etaValue then
      (__old_a,__old_eta) := aValue,etaValue;
      ode_sol('parameters'=[a=aValue,eta=etaValue]);
   end if;
   res:=rhs(ode_sol(tValue)[2]);
   #evalf[8](res);
end proc:
> 

Be11Num(a,eta,t); # returns unevaluated, good

Be11Num(a, eta, t)

> 

ode_sol(parameters=[a=1,eta=0.1]);
ode_sol(2.0)[2], Be11Num(1,0.1,2.0); # these should agree

[a = 1., eta = .1]

Be11(t) = HFloat(1.2342315801626595), HFloat(1.2342315801626595)

> 

ode_sol(parameters=[a=1.4,eta=0.3]);
ode_sol(2.0)[2], Be11Num(1.4,0.3,2.0); # these should agree

[a = 1.4, eta = .3]

Be11(t) = HFloat(1.1815627486067706), HFloat(1.1815627486067706)

> 

tensTemp1:=2*C1*tensB+4*C2*(Trace(tensB)-3)*tensB+6*C3*(Trace(tensB)-3)*tensB-2*C4*MatrixInverse(tensB):
tensTemp1Dev:=tensTemp1-(1/3)*Trace(tensTemp1)*Matrix(3,shape=identity):
tensBe:=Matrix([[Be11Num(a,eta,t),0,0],[0,1/sqrt(Be11Num(a,eta,t)),0],[0,0,1/sqrt(Be11Num(a,eta,t))]]):
tensTemp2Dev:=2*a*tensBe-(1/3)*Trace(2*a*tensBe)*~Matrix(3,shape=identity):
sigma:=tensTemp1Dev+tensTemp2Dev-p*~Matrix(3,shape=identity):
p:=sigma[2,2]+p:
sigma:=tensTemp1Dev+tensTemp2Dev-p*~Matrix(3,shape=identity):
PK:=sigma.Transpose(MatrixInverse(tensF)):
PK:=PK[1,1]: #expression  to fit
> 

vectorT:=[seq(i,i=0..13/3,0.01)]:
vectorT:=Vector(vectorT):
> 

DS2:=CodeTools:-Usage(
        DirectSearch:-DataFit(PK,vectorT,PK1Expe,t,
                              [ a=0.1 .. 10.0, eta=0.1 .. 10.2
                                ,C1=0..1, C2=0..10, C3=-6..6, C4=0..1
                              ])
);

memory used=7.16GiB, alloc change=6.00MiB, cpu time=110.76s, real time=105.76s, gc time=10.59s

[HFloat(5.798836840043786e-7), [C1 = HFloat(0.07091614374135646), C2 = HFloat(0.19080181427660822), C3 = HFloat(-0.13066108619768044), C4 = HFloat(2.090135897538898e-5), a = HFloat(0.6139583977913079), eta = HFloat(0.5024357747664052)], 615]

> 

plots:-display(
   plot(<vectorT|PK1Expe>, style=point, color=blue),
   plot(eval(PK,DS2[2]), t=min(vectorT)..max(vectorT),
   size=[600,300]));

> 

 

Download donnees_cycles_sin_DS.mw

First 109 110 111 112 113 114 115 Last Page 111 of 336