acer

32333 Reputation

29 Badges

19 years, 323 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

That is the return value of the call to Tabulate.

You can simply use a statement-terminating colon to suppress it.

Why not write a shell script (text file, with suitable hash bang first line), and ensure that sets its own environment adequately? You could pass arguments along to that, via ssystem.

It's not clear why all you tried before lead nowhere. Why are .cshrc and tcsh and python not referenced with absolute paths?

Surely you should instead refer to shell executable by explicit full location (whether it be to a new shell script file, or to tcsh, or what else that is not a 'sh' builtin).

And how you form strings of commands (possible containing escaped double quotes, or other sttings) to pass along  can matter. It's not clear precisely what you have tried.

 

You asked at least one other question in the last two weeks about patmatch, where the convincing answer was that indets is the better tool for the job at hand.

In at least one exchange you seemed to resist that idea.

And now you seem to have completely forgotten about indets altogether. (That seems weird, to me.)

Please consider the following: indets is one of the most basic tools in Maple, being both fast and flexible. As a result it is very commonly used within stock Library commands.

In contrast patmatch can be slow and quite difficult to use flexibly and robustly. It is very rarely used in stock Library commands.

Your current example is handled by a simple call to indets. I encourage you to figure out the exact call yourself, as doing so would help your understanding and could greatly aid your future Maple development work.

note: the great advantage of indets is that it works with Maple types. The type system is very highly developed in Maple and works closely with the underlying structure of expressions. It's understandable that a Mathematica user would think first to code with patterns, since Mma's pattern-matching is highly developed and works more closely in conjunction with its expression structures. Since Maple and Mathematica have different expression strucures (and supporting tools) then the best mechanisms for matching substructures are also different.

Some years ago I answered a similar question, using a hand-rolled procedure I named &// which would apply its second argument to its first.

Eg,

   foo &// simplify;

would execute simplify(foo). This syntax works even if foo is some involved function call -- where foo gets applied to the result of that.

I only have a cell phone right now, and cannot easily find that old post. I'm sure someone else could write it quickly.

I expect it went something like this (possibly with a few defensive checks)

`&//` := proc(a,b)
             b(a):
end proc;

So then it could get used as follows,

cos(x)^2+sin(x)^2 &// simplify;

Sometimes (but not always), conversion to exp form can help. But the effort to simplify may also vary by example.

And the Roots command can sometime be useful for results within a finite range.

It just so happens that convert to exp form helps produce a simpler set of parametrized formula for the very first example. But I have seen examples where conversion to exp form prevent solve from finding the proper full general solution.

Below I include both the original and the revised example. The techniques that can be practically applied are slightly different.

restart;

kernelopts(version);

`Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874`

eq := sin(9*x-(1/3)*Pi) = sin(7*x-(1/3)*Pi):

S2 := [solve(convert({eq},exp), x, allsolutions, explicit=true)]:

SS2 := {map(s->eval(x,s),simplify(S2))[]};

{-Pi*_Z1, -(19/48)*Pi-Pi*_Z3, -(13/48)*Pi-Pi*_Z2, -(7/48)*Pi-Pi*_Z3, -(1/48)*Pi-Pi*_Z2, (5/48)*Pi-Pi*_Z3, (11/48)*Pi-Pi*_Z2, (17/48)*Pi-Pi*_Z3, (23/48)*Pi-Pi*_Z2}

pts := {Student:-Calculus1:-Roots(eq, x=-2*Pi..2*Pi, maxsols=100)[]};

{0, Pi, -2*Pi, -Pi, 2*Pi, -(91/48)*Pi, -(85/48)*Pi, -(79/48)*Pi, -(73/48)*Pi, -(67/48)*Pi, -(61/48)*Pi, -(55/48)*Pi, -(49/48)*Pi, -(43/48)*Pi, -(37/48)*Pi, -(31/48)*Pi, -(25/48)*Pi, -(19/48)*Pi, -(13/48)*Pi, -(7/48)*Pi, -(1/48)*Pi, (5/48)*Pi, (11/48)*Pi, (17/48)*Pi, (23/48)*Pi, (29/48)*Pi, (35/48)*Pi, (41/48)*Pi, (47/48)*Pi, (53/48)*Pi, (59/48)*Pi, (65/48)*Pi, (71/48)*Pi, (77/48)*Pi, (83/48)*Pi, (89/48)*Pi, (95/48)*Pi}

plots:-display([
        plot((rhs-lhs)(eq),x=-2*Pi..2*Pi,color=blue),
        plots:-pointplot(map(p->[p,eval((rhs-lhs)(eq),x=p)],pts),
                         color=red,symbol=circle,symbolsize=20)],
        size=[700,200]
       );

iSS2:=indets(SS2,name) minus {Pi}:
pts2:=select(r->is(r>=-2*Pi and r<=2*Pi),
       `union`(seq(seq(seq(eval(SS2,[iSS2[1]=i,iSS2[2]=j,iSS2[3]=k]),i=-2..2),j=-2..2),k=-2..2)));

{0, Pi, -2*Pi, -Pi, 2*Pi, -(91/48)*Pi, -(85/48)*Pi, -(79/48)*Pi, -(73/48)*Pi, -(67/48)*Pi, -(61/48)*Pi, -(55/48)*Pi, -(49/48)*Pi, -(43/48)*Pi, -(37/48)*Pi, -(31/48)*Pi, -(25/48)*Pi, -(19/48)*Pi, -(13/48)*Pi, -(7/48)*Pi, -(1/48)*Pi, (5/48)*Pi, (11/48)*Pi, (17/48)*Pi, (23/48)*Pi, (29/48)*Pi, (35/48)*Pi, (41/48)*Pi, (47/48)*Pi, (53/48)*Pi, (59/48)*Pi, (65/48)*Pi, (71/48)*Pi, (77/48)*Pi, (83/48)*Pi, (89/48)*Pi, (95/48)*Pi}

pts2 minus pts, pts minus pts2;

{}, {}

restart;

eq := sin(9*x-(1/3)*Pi) = sin(5*x-(1/6)*Pi);

-cos(9*x+(1/6)*Pi) = -cos(5*x+(1/3)*Pi)

S2 := [solve(convert({eq},exp), x, allsolutions, explicit=true)]:

SS2 := {map(s->eval(x,s),convert(simplify(evalc(simplify(S2))),tan))[]};

{-(13/28)*Pi+Pi*_Z2, -(11/24)*Pi+Pi*_Z1, -(9/28)*Pi+Pi*_Z2, -(5/28)*Pi+Pi*_Z2, -(1/28)*Pi+Pi*_Z2, (1/4)*Pi+Pi*_Z2, (1/24)*Pi+Pi*_Z1, (3/28)*Pi+Pi*_Z2, (11/28)*Pi+Pi*_Z2}

pts := identify({Student:-Calculus1:-Roots(eq, x=-2*Pi..2*Pi, maxsols=100, numeric)[]});

{-(53/28)*Pi, -(47/24)*Pi, -(45/28)*Pi, -(41/28)*Pi, -(37/28)*Pi, -(35/24)*Pi, -(33/28)*Pi, -(29/28)*Pi, -(25/28)*Pi, -(23/24)*Pi, -(17/28)*Pi, -(13/28)*Pi, -(11/24)*Pi, -(9/28)*Pi, -(7/4)*Pi, -(5/28)*Pi, -(3/4)*Pi, -(1/28)*Pi, (1/4)*Pi, (1/24)*Pi, (3/28)*Pi, (5/4)*Pi, (11/28)*Pi, (13/24)*Pi, (15/28)*Pi, (19/28)*Pi, (23/28)*Pi, (25/24)*Pi, (27/28)*Pi, (31/28)*Pi, (37/24)*Pi, (39/28)*Pi, (43/28)*Pi, (47/28)*Pi, (51/28)*Pi, (55/28)*Pi}

plots:-display([
        plot((rhs-lhs)(eq),x=-2*Pi..2*Pi,color=blue),
        plots:-pointplot(map(p->[p,eval((rhs-lhs)(eq),x=p)],pts),
                         color=red,symbol=circle,symbolsize=20)],
        size=[700,200]
       );

iSS2:=indets(SS2,name) minus {Pi}:
pts2:=select(r->is(r>=-2*Pi and r<=2*Pi),
       `union`(seq(seq(seq(eval(SS2,[iSS2[1]=i,iSS2[2]=j]),i=-2..2),j=-2..2))));

{-(53/28)*Pi, -(47/24)*Pi, -(45/28)*Pi, -(41/28)*Pi, -(37/28)*Pi, -(35/24)*Pi, -(33/28)*Pi, -(29/28)*Pi, -(25/28)*Pi, -(23/24)*Pi, -(17/28)*Pi, -(13/28)*Pi, -(11/24)*Pi, -(9/28)*Pi, -(7/4)*Pi, -(5/28)*Pi, -(3/4)*Pi, -(1/28)*Pi, (1/4)*Pi, (1/24)*Pi, (3/28)*Pi, (5/4)*Pi, (11/28)*Pi, (13/24)*Pi, (15/28)*Pi, (19/28)*Pi, (23/28)*Pi, (25/24)*Pi, (27/28)*Pi, (31/28)*Pi, (37/24)*Pi, (39/28)*Pi, (43/28)*Pi, (47/28)*Pi, (51/28)*Pi, (55/28)*Pi}

pts2 minus pts, pts minus pts2;

{}, {}

 

Download trig_sols2.mw

(I wrote this before Mariusz's followup for the second example, which also does conversion to exp form.)

Here's one way. (No doubt it could be slicker.)

restart;

blockify:=proc(M::Matrix,m::posint)
    Matrix(convert(map(e->LinearAlgebra:-IdentityMatrix(m)*e, M),listlist));
end proc:

A:= Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 2});

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

blockify(A, 2);

Matrix(4, 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) = 2, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 2})

B := Matrix([[-7,0],[0,k]]);

Matrix(2, 2, {(1, 1) = -7, (1, 2) = 0, (2, 1) = 0, (2, 2) = k})

blockify(B, 3);

Matrix(6, 6, {(1, 1) = -7, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = -7, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -7, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = k, (4, 5) = 0, (4, 6) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = k, (5, 6) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = k})

C := LinearAlgebra:-RandomMatrix(3, generator=-1..1);

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

blockify(C, 2);

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

 

Download blockify.mw

First I downloaded that problematic Dtm.mw file, and then I ran the following command using the location to which I saved it.

The generated string of 1D code is not perfect (it has some escaped double-quotes). But it visually shows S*(k+1) as the left hand side of an invalid assignment. And if I manually remove the escaped double-quotes (escaped special characters) and try to run it then I receive an error message about that as an invalid assignment.

restart;

Worksheet:-WorksheetToMapleText(cat(kernelopts(homedir),"/mapleprimes/Dtm.mw"));

"restart: _local(gamma); m := 10; A := 10; delta := .112; rho := .23; mu := .31; beta := 1.4; alpha := 2.1; gamma := 1.02; q := 2.3; b1 := 50; b2 := 10; b3 := 5; b4 := 20; S(0) := b1; B(0) := b2; V(0) := b3; R(0) := b4; "   for k from 0 to N do "S*(k+1) := 1/(k+1)*[A*delta*k-(rho+`&micro;`)*S*k-beta*sum(S*m*B*(k-m),k = 0 .. m)]; B*(k+1) := 1/(k+1)*[(`&micro;`+alpha+gamma)*B*k-beta*sum(S*m*B*(k-m),k = 0 .. m)]; V*(k+1) := 1/(k+1)*[rho*S*k-(1-q)*S*k-`&micro;`*V(k)]; R*(k+1) := 1/(k+1)*[-`&micro;`*R*k+B*gamma*k]; "  end do;";
"

 

Download Dtm_text.mw

[edited] It's understandable, perhaps, that this utility might baulk at complete conversion of input that has a syntax error. I think that it's obvious that code containing a syntax error is in general not translatable. Even though the incorrect syntax for the given example does have a rather obvious and natural 1D analogue.

You might concoct your own procedure to instantiate a given RV at specific values for its parameters.

For example (and you could make this more bullet proof with some additional, defensive checks),

restart;

instantiate:=proc(r::name,L::list(name=anything))
  uses ST=Statistics;
  local d;
  d := [attributes(X)][3]; # you could make this defensive
  ST:-RandomVariable(d:-ParentName(op(eval(d:-Parameters,L))));
end proc:

 

X := Statistics:-RandomVariable(Normal(mu, sigma));

_R

Y := instantiate(X, [mu=1]);

_R0

Statistics:-Mean(Y);

1

Statistics:-PDF(Y,t);

(1/2)*2^(1/2)*exp(-(1/2)*(t-1)^2/sigma^2)/(Pi^(1/2)*sigma)

 

Download RV_instantiate.mw

This also works for deeper nesting.

A1 := {1,2,{5,6},3,{13,{{21,27},16}},{8,9},{{11,12},4}};

        A1 := {1, 2, 3, {4, {11, 12}}, {5, 6},
               {8, 9},{13, {16, {21, 27}}}}

{subsindets(A1,set,op)};

        {1, 2, 3, 4, 5, 6, 8, 9, 11, 12, 13, 16, 21, 27}
restart;

X := Statistics:-RandomVariable(Normal(0, 1)):
pp:=[attributes(X)][3]:
pp:-ParentName;
                       NormalDistribution
pp:-Parameters;
                             [0, 1]

Y := Statistics:-RandomVariable(Normal(5, 17)):
pp:=[attributes(Y)][3]:
pp:-ParentName;
                       NormalDistribution
pp:-Parameters;
                            [5, 17]

Z := Statistics:-RandomVariable(Exponential(3)):
pp:=[attributes(Z)][3]:
pp:-ParentName;
                    ExponentialDistribution
pp:-Parameters;
                              [3]

You can easily construct your own float[8] Array to denote the color, using either HUE or RGB specification.

restart;

x1:=[seq(x,x=1..2,0.02)]:
y1:=[seq(y^3,y=1..2,0.02)]:

N:=nops(x1);
z1:=[seq(x1[i]/y1[i],i=1..N)]:

51

#
# HUE requires one float per point, to denote the color.
#
plots:-pointplot([x1,y1],
                 gridlines=false,
                 symbol=solidcircle, symbolsize=15,
                 color=COLOR(HUE,Array(1..N,z1,datatype=float[8])));

#
# RGB requires three floats per point, to denote the color.
#
C:=Array(1..51,1..3,datatype=float[8]):
#
# Alter only blue layer here.
#
C[..,3]:=Array(1..51,z1,datatype=float[8]):

plots:-pointplot([x1,y1],
                 gridlines=false,
                 symbol=solidcircle, symbolsize=15,
                 color=COLOR(RGB,C));

#
# Alter blue and green layers here (differently).
#
C[..,3]:=map(c->1-(c-min(z1))/(max(z1)-min(z1)), Array(1..51,z1,datatype=float[8])):
C[..,2]:=map(c->(c-min(z1))/(max(z1)-min(z1)), Array(1..51,z1,datatype=float[8])):
plots:-pointplot([x1,y1],
                 gridlines=false,
                 symbol=solidcircle, symbolsize=15,
                 color=COLOR(RGB,C));

 

Download color_pointplot.mw

On MS-Windows I add the following undocumented option (following a blank space) to the launch command in the Properties of the desktop icon,

  -standalone

The main effect, apparently, is that the icon subsequently launches wholly separate java virtual machines.

If one of these GUI instances crashes (while printing a huge result, rotating a huge plot, waiting on a crashed Maple kernel, etc) then the other GUI instance remains responsive and its worksheets can still be saved properly and promptly, etc.

Since I quite often run multiple GUI sessions, often with examples that can crash the GUI, this has been very useful in saving my other work in progress.

My machine has adequate RAM, and I don't care that separate GUI JVMs can take up 300MB footprint each while running. It is much more valuable to me that I don't lose work in progress.

(Sorry, I did not notice this Question before now.)

 

One way to disable the discontinuity check done by evalf(Int(...)) is to use operator form for the integrand. And fsolve prefers everything in operator form when it detects such (even if not really necessary for a numeric scalar like 0.5 on the rhs).

restart;
fe := Int(x -> sqrt(1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)
               -5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2), 0 .. al) = 0.5:

fsolve(unapply((lhs-rhs)(fe), al), 0..1);

                          0.1550088024

That works even when the original int call had already been formed for the given example: We can recast the int call with integrand in expression form to an Int call with integrand in operator form.

restart;
fe := int(sqrt(1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)
          -5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2), x=0 .. al) = 0.5:

F := subsindets(fe, specfunc({int,Int}),
                u->Int(unapply(op(1,u),lhs(op(2,u))),rhs(op(2,u)),op(3..,u))):

fsolve(unapply((lhs-rhs)(F), al), 0..1);

                          0.1550088024

Now some comments about Maple's numeric integration.

It is long-standing and unfortunate behaviour that evalf(Int(...)) (ie, numeric integration) offers no all-round easy option to disable its examination of the integrand for discontinuities.

The important 1D numeric quadrature schemes utilized have error estimates that rely on certain degrees of smoothness, and thus they should often split the domain in the presence of a finite number of discontinuities. But there are lots of examples where the check for discontinuity consumes enormous resources, and where that check is not necessary or wanted (because there are no such points).

The typical unfortunate scenario occurs when a symbolic discontinutiy check on an expression containing float coefficients causes a symbolic expansion which runs amok when the floats get converted to exact rationals with large numerator or denominator.

Sometimes the check can be avoided by specifiying the method=_d01ajc for evalf(Int(...)) example. But that method requires Digits<=15, and when the calling environment has default Digits=10 then fsolve may try and resolve the question of convergence (or a root candidate) at yet higher working precision.

While it would indeed be useful if fsolve were to accept an accuracy tolerance option, the runaway crash in this particular example is due (IMO) to a bug in evalf(Int(...)) and not in fsolve.

In my opinion it is a very longstanding and serious bug that evalf(Int(...)) offers no option to forcibly disable its discontinuity test. I have answered dozens of Mapleprimes questions over the years where evalf(Int(...)) showed serious performance problems for users, where it was not originally obvious that the discontinuity check was the culprit.

I have seen several examples that don't involve fsolve at all, where the evalf(Int(...)) misbehaviour causes unnecessary very slow performance or even a crash. Some of these are plotting examples. Some are simply single instances of a numeric integration. The central problem is not in fsolve.

Consider this example, which runs amok (don't try it if you don't want a crash):

restart;
Q := Int(sqrt(1+(-5.557990765*sin(5.557990765*x)-7.3*cos(5.557990765*x)
         -5.6*sinh(5.557990765*x)+7.3*cosh(5.557990765*x))^2), x=0 .. al):
evalf[25](eval(Q,al=0.3));

The way that I disable the discontinuity check at the start of this Answer is to use operator form for the integrand. That can become very tricky when the problem involves nested integrals or integrands that depend on (scoped) parameters, etc. It can become so tricky that it really cannot be considered a practical alternative to a (currently missing) dedicated option to prevent the discontinuity check.

Here is the first plot obtained in the worksheet below, but without the blue plane. It has a smooth lower edge and no missing portion (for larger y values).

An upper bound y=20 was supplied. For the first approach below you could increase that (and the `view`, accordingly). An upper bound for the z was also supplied, in the `view` option, which you could also adjust to taste.

restart;

solve({x^3*y*z=1/2,x<y , y<z}, [z], real, parametric, explicit);

piecewise(x < 0 and x < y and y < 0, [[z = 1/(2*y*x^3)]], `and`(`and`(`and`(0 < x, x < RootOf(2*_Z^5-1, .8125 .. .875)), x < y), y < sqrt(2)*sqrt(x^3)/(2*x^3)), [[z = 1/(2*y*x^3)]], [])

#
# You can exclude the blue plane if you desire.
#
plots:-display(
  plot3d(1/(2*y*x^3), x=0 .. RootOf(2*_Z^5-1,0.8125..0.875),
                      y=x..sqrt(2)*sqrt(x^3)/(2*x^3),
         color=green),
  plot3d(y, x=0..1, y=0..20,
         color=blue),
  view=[0..1,0..20,0..100]
);

#
# One problem with this approach is that a portion of the
# green surface is missing (and the missing portion is
# reduced only as the grid resolution is made very fine).
#
# Another problem is that the blue plane is needed to
# delimit the correct green surface's lower boundary.
#
f:=x^3*y*z:
S:=solve(f=1/2, z);
plot3d([S, y], x=0..y, y=0..20,
       color=[green,blue],grid=[200,200],
       view=[0..1,0..20,0..100]);

(1/2)/(y*x^3)

#
# It is not easy to produce a smooth lower edge for the
# valid part of the green surface.
#
# This also increases the missing portion of the green
# surface for larger y values.
#
plot3d([piecewise(S>y,S,undefined), y], x=0..y, y=0..20,
       color=[green,blue],grid=[200,200],
       view=[0..1,0..20,0..100]);

 

Download plot3d_inequality.mw

Here are a few ways, each involving just a few calls to LinearAlgebra:-GenerateMatrix.

I didn't keep exactly the same names as you have for Vector X, etc.

I include the situation where the final LHS is left in the form Matrix.diff(X) , as opposed to being fully explicit.

(But I would not be surprised if there were utilities to help in DEtools or PDEtools.)

[edit] I have included a way to print the Matrix equations (using either X.DX, or DX explicitly on the LHS) with inert multiplication. In case you were after such a visual display. For the case of the explicit LHS I've separated the inverse of X using its determinant factor and its adjoint.

restart;

eq1:=(L1/n12 + n12*L2)*diff(i2(t), t) + L1*diff(i3(t), t)/n13 = -v1(t) - R1*i3(t)/n13 - R1*i2(t)/n12 + n12*v2(t) - n12*R2*i2(t);

eq2:=L1*diff(i2(t), t)/n12 + (L1/n13 + n13*L3)*diff(i3(t), t) = -v1(t) - R1*i3(t)/n13 - R1*i2(t)/n12 + n13*v3(t) - n13*R3*i3(t);

eq3:=-L2*diff(i2(t), t) + L1*diff(i1(t), t)/n12 = -v2(t) + R2*i2(t) + v1(t)/n12 - R1*i1(t)/n12;

(L1/n12+n12*L2)*(diff(i2(t), t))+L1*(diff(i3(t), t))/n13 = -v1(t)-R1*i3(t)/n13-R1*i2(t)/n12+n12*v2(t)-n12*R2*i2(t)

L1*(diff(i2(t), t))/n12+(L1/n13+n13*L3)*(diff(i3(t), t)) = -v1(t)-R1*i3(t)/n13-R1*i2(t)/n12+n13*v3(t)-n13*R3*i3(t)

-L2*(diff(i2(t), t))+L1*(diff(i1(t), t))/n12 = -v2(t)+R2*i2(t)+v1(t)/n12-R1*i1(t)/n12

with(LinearAlgebra):

DT:=[diff(i1(t),t),diff(i2(t),t),diff(i3(t),t)]:
DX:=Vector(DT):

iT:=[i1(t),i2(t),i3(t)]:
Vi:=Vector(iT):

vT:=[v1(t),v2(t),v3(t)]:
Vv:=Vector(vT):

X:=GenerateMatrix(map(lhs,[eq1,eq2,eq3]),DT)[1];
altA:=GenerateMatrix(map[2](select,has,expand(map(rhs,[eq1,eq2,eq3])),iT),iT)[1];
altB:=GenerateMatrix(map[2](select,has,expand(map(rhs,[eq1,eq2,eq3])),vT),vT)[1];

Matrix(3, 3, {(1, 1) = 0, (1, 2) = L1/n12+n12*L2, (1, 3) = L1/n13, (2, 1) = 0, (2, 2) = L1/n12, (2, 3) = L1/n13+n13*L3, (3, 1) = L1/n12, (3, 2) = -L2, (3, 3) = 0})

Matrix(3, 3, {(1, 1) = 0, (1, 2) = -R1/n12-n12*R2, (1, 3) = -R1/n13, (2, 1) = 0, (2, 2) = -R1/n12, (2, 3) = -R1/n13-n13*R3, (3, 1) = -R1/n12, (3, 2) = R2, (3, 3) = 0})

Matrix(%id = 18446884362352865638)

X.DX = map(expand, altA.Vi + altB.Vv );

(Vector(3, {(1) = (L1/n12+n12*L2)*(diff(i2(t), t))+L1*(diff(i3(t), t))/n13, (2) = L1*(diff(i2(t), t))/n12+(L1/n13+n13*L3)*(diff(i3(t), t)), (3) = -L2*(diff(i2(t), t))+L1*(diff(i1(t), t))/n12})) = (Vector(3, {(1) = -v1(t)-R1*i3(t)/n13-R1*i2(t)/n12+n12*v2(t)-n12*R2*i2(t), (2) = -v1(t)-R1*i3(t)/n13-R1*i2(t)/n12+n13*v3(t)-n13*R3*i3(t), (3) = -v2(t)+R2*i2(t)+v1(t)/n12-R1*i1(t)/n12}))

T:=map(rhs,solve([eq1,eq2,eq3],DT)[1]):

A:=GenerateMatrix(map[2](select,has,expand(T),iT),iT)[1]:

B:=GenerateMatrix(map[2](select,has,expand(T),vT),vT)[1]:

simplify(DX = A.Vi + B.Vv);

(Vector(3, {(1) = diff(i1(t), t), (2) = diff(i2(t), t), (3) = diff(i3(t), t)})) = (Vector(3, {(1) = (-(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2)*R1*i1(t)+L3*n12*n13^2*(L1*R2-L2*R1)*i2(t)+L2*n12^2*n13*(L1*R3-L3*R1)*i3(t)+L1*((L2*n12^2+L3*n13^2)*v1(t)-v2(t)*L3*n12*n13^2-v3(t)*L2*n12^2*n13))/(L1*(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2)), (2) = ((-R2*(L3*n13^2+L1)*n12^2-L3*R1*n13^2)*i2(t)+(n13*(L1*R3-L3*R1)*i3(t)+n12*(L3*n13^2+L1)*v2(t)-v1(t)*L3*n13^2-v3(t)*L1*n13)*n12)/(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2), (3) = ((-R3*(L2*n12^2+L1)*n13^2-L2*R1*n12^2)*i3(t)+(n12*(L1*R2-L2*R1)*i2(t)+n13*(L2*n12^2+L1)*v3(t)-v1(t)*L2*n12^2-v2(t)*L1*n12)*n13)/(L3*(L2*n12^2+L1)*n13^2+L1*L2*n12^2)}))

DX = simplify(LinearSolve(X, map(expand, altA.Vi + altB.Vv )));

(Vector(3, {(1) = diff(i1(t), t), (2) = diff(i2(t), t), (3) = diff(i3(t), t)})) = (Vector(3, {(1) = (-(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2)*R1*i1(t)+L3*n12*n13^2*(L1*R2-L2*R1)*i2(t)+L2*n12^2*n13*(L1*R3-L3*R1)*i3(t)+L1*((L2*n12^2+L3*n13^2)*v1(t)-v2(t)*L3*n12*n13^2-v3(t)*L2*n12^2*n13))/(L1*(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2)), (2) = ((-R2*(L3*n13^2+L1)*n12^2-L3*R1*n13^2)*i2(t)+(n13*(L1*R3-L3*R1)*i3(t)+n12*(L3*n13^2+L1)*v2(t)-v1(t)*L3*n13^2-v3(t)*L1*n13)*n12)/(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2), (3) = ((-R3*(L2*n12^2+L1)*n13^2-L2*R1*n12^2)*i3(t)+(n12*(L1*R2-L2*R1)*i2(t)+n13*(L2*n12^2+L1)*v3(t)-v1(t)*L2*n12^2-v2(t)*L1*n12)*n13)/(L3*(L2*n12^2+L1)*n13^2+L1*L2*n12^2)}))

Xinv:=X^(-1):
fullA:=simplify(Xinv.altA):
fullB:=simplify(Xinv.altB):

DX = simplify(fullA.Vi + fullB.Vv);

(Vector(3, {(1) = diff(i1(t), t), (2) = diff(i2(t), t), (3) = diff(i3(t), t)})) = (Vector(3, {(1) = (-(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2)*R1*i1(t)+L3*n12*n13^2*(L1*R2-L2*R1)*i2(t)+L2*n12^2*n13*(L1*R3-L3*R1)*i3(t)+L1*((L2*n12^2+L3*n13^2)*v1(t)-v2(t)*L3*n12*n13^2-v3(t)*L2*n12^2*n13))/(L1*(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2)), (2) = ((-R2*(L3*n13^2+L1)*n12^2-L3*R1*n13^2)*i2(t)+(n13*(L1*R3-L3*R1)*i3(t)+n12*(L3*n13^2+L1)*v2(t)-v1(t)*L3*n13^2-v3(t)*L1*n13)*n12)/(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2), (3) = ((-R3*(L2*n12^2+L1)*n13^2-L2*R1*n12^2)*i3(t)+(n12*(L1*R2-L2*R1)*i2(t)+n13*(L2*n12^2+L1)*v3(t)-v1(t)*L2*n12^2-v2(t)*L1*n12)*n13)/(L3*(L2*n12^2+L1)*n13^2+L1*L2*n12^2)}))

simplify(A - fullA), simplify(B - fullB);

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}), 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})

X %. DX = altA %. Vi + altB %. Vv;

`%.`(Matrix(3, 3, {(1, 1) = 0, (1, 2) = L1/n12+n12*L2, (1, 3) = L1/n13, (2, 1) = 0, (2, 2) = L1/n12, (2, 3) = L1/n13+n13*L3, (3, 1) = L1/n12, (3, 2) = -L2, (3, 3) = 0}), Vector(3, {(1) = diff(i1(t), t), (2) = diff(i2(t), t), (3) = diff(i3(t), t)})) = `%.`(Matrix(3, 3, {(1, 1) = 0, (1, 2) = -R1/n12-n12*R2, (1, 3) = -R1/n13, (2, 1) = 0, (2, 2) = -R1/n12, (2, 3) = -R1/n13-n13*R3, (3, 1) = -R1/n12, (3, 2) = R2, (3, 3) = 0}), Vector(3, {(1) = i1(t), (2) = i2(t), (3) = i3(t)}))+`%.`(Matrix(3, 3, {(1, 1) = -1, (1, 2) = n12, (1, 3) = 0, (2, 1) = -1, (2, 2) = 0, (2, 3) = n13, (3, 1) = 1/n12, (3, 2) = -1, (3, 3) = 0}), Vector(3, {(1) = v1(t), (2) = v2(t), (3) = v3(t)}))

DX = ((1/Determinant(X)) %* Adjoint(X)) %. (altA %. Vi + altB %. Vv);

(Vector(3, {(1) = diff(i1(t), t), (2) = diff(i2(t), t), (3) = diff(i3(t), t)})) = `%.`(`%*`(n12^2*n13/(L1*(L2*L3*n12^2*n13^2+L1*L2*n12^2+L1*L3*n13^2)), Matrix(3, 3, {(1, 1) = L2*(L3*n13^2+L1)/n13, (1, 2) = -L1*L2/n13, (1, 3) = (L2*L3*n12^2*n13^2+L1*L2*n12^2+L1*L3*n13^2)/(n12*n13), (2, 1) = (L3*n13^2+L1)*L1/(n12*n13), (2, 2) = -L1^2/(n12*n13), (2, 3) = 0, (3, 1) = -L1^2/n12^2, (3, 2) = (L2*n12^2+L1)*L1/n12^2, (3, 3) = 0})), `%.`(Matrix(3, 3, {(1, 1) = 0, (1, 2) = -R1/n12-n12*R2, (1, 3) = -R1/n13, (2, 1) = 0, (2, 2) = -R1/n12, (2, 3) = -R1/n13-n13*R3, (3, 1) = -R1/n12, (3, 2) = R2, (3, 3) = 0}), Vector(3, {(1) = i1(t), (2) = i2(t), (3) = i3(t)}))+`%.`(Matrix(3, 3, {(1, 1) = -1, (1, 2) = n12, (1, 3) = 0, (2, 1) = -1, (2, 2) = 0, (2, 3) = n13, (3, 1) = 1/n12, (3, 2) = -1, (3, 3) = 0}), Vector(3, {(1) = v1(t), (2) = v2(t), (3) = v3(t)})))

simplify(value(%));

(Vector(3, {(1) = diff(i1(t), t), (2) = diff(i2(t), t), (3) = diff(i3(t), t)})) = (Vector(3, {(1) = (-(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2)*R1*i1(t)+L3*n12*n13^2*(L1*R2-L2*R1)*i2(t)+L2*n12^2*n13*(L1*R3-L3*R1)*i3(t)+L1*((L2*n12^2+L3*n13^2)*v1(t)-v2(t)*L3*n12*n13^2-v3(t)*L2*n12^2*n13))/(L1*(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2)), (2) = ((-R2*(L3*n13^2+L1)*n12^2-L3*R1*n13^2)*i2(t)+(n13*(L1*R3-L3*R1)*i3(t)+n12*(L3*n13^2+L1)*v2(t)-v1(t)*L3*n13^2-v3(t)*L1*n13)*n12)/(L2*(L3*n13^2+L1)*n12^2+L1*L3*n13^2), (3) = ((-R3*(L2*n12^2+L1)*n13^2-L2*R1*n12^2)*i3(t)+(n12*(L1*R2-L2*R1)*i2(t)+n13*(L2*n12^2+L1)*v3(t)-v1(t)*L2*n12^2-v2(t)*L1*n12)*n13)/(L3*(L2*n12^2+L1)*n13^2+L1*L2*n12^2)}))

 

Download various2.mw

First 148 149 150 151 152 153 154 Last Page 150 of 336