acer

32333 Reputation

29 Badges

19 years, 320 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

OpProp is a procedure. Procedures have last-name-eval, so when you invoke just,
   lastexception;
you just get the expression sequence containing the literal name OrProp. But when you do
   print(lastexception);
you are requesting more (depending on the value of interface(verboseproc) at the time).

You could also compare the behaviour for,
   lastexception[1];
versus,
   eval(lastexception[1]);
and,
   print(OrProp);
It is unusual, that the last example immediately above also shows the issue.

In the CommandLine interface your code results in the body of OrProp being printed to the terminal. It's not that long.

I don't know yet why the Java GUI thinks that it has to elide altogether the printing of OrProp as a procedure. But I notice a discrepency between the results from,
   printf("%P",eval(OrProp));
and that of,
   showstat(OrProp);
around the code line numbers 52 and 53. It might be related.

But it's stranger because the behavior you've mentioned doesn't seem to depend on the setting of interface(verboseproc) .

I know that you are interested in using timelimit. You could try something like this, instead,

   print(sprintf("%a",lastexception[1]),
            StringTools:-FormatMessage(lastexception[2..-1]));

I mentioned that last bit as I suspect you might not be fully satisfied with the results from your suggestion,
   StringTools:-FormatMessage(lastexception)

Other useful tools include stopat, trace, and showstat.

For example, using trace (after using stopat and showstat to discover the computation path),

restart;

infolevel[dsolve]:=6;

6

ode:=x^2*diff(diff(y(x),x),x)+x*diff(y(x),x)+(x^2-1/4)*y(x)=0;

x^2*(diff(diff(y(x), x), x))+x*(diff(y(x), x))+(x^2-1/4)*y(x) = 0

#trace(`ODEtools/odsolve`):
#trace(`odsolve/2nd_order`);
#trace(`odsolve/2nd_order/class`);
#trace(`odsolve/2nd_order/linear_homogeneous`);

##trace(`odsolve/2nd_order/quadrature`);
##trace(`odsolve/2nd_order/const_coeffs`);
##trace(`odsolve/2nd_order/Euler`);
##trace(`odsolve/2nd_order/linear_1`);
##trace(`odsolve/2nd_order/linear/missing_y`);

#trace(`odsolve/2nd_order/Kovacic`);
#trace(`odeadv/2nd_order/Kovacic`);
#trace(`DEtools/kovacicsols`);
trace(`DEtools/kovacicsols/order2`);

`DEtools/kovacicsols/order2`

dsolve(ode);

Methods for second order ODEs:

--- Trying classification methods ---

trying a quadrature

checking if the LODE has constant coefficients

checking if the LODE is of Euler type

trying a symmetry of the form [xi=0, eta=F(x)]

   testing BRANCH 1 ->

   testing BRANCH 2 ->

   testing BRANCH 3 ->

   testing BRANCH 4 ->

   testing BRANCH 5 ->

   testing BRANCH 6 ->

checking if the LODE is missing 'y'

-> Trying a Liouvillian solution using Kovacic's algorithm

{--> enter \`DEtools/kovacicsols/order2\`, args = [4*x^2-1, 4*x, 4*x^2], x

{--> enter \`DEtools/kovacicsols/order2\`, args = [4*x^2-1, 4*x, 4*x^2], x, {}

{}, {}, {}, {}

[1, 2, [4, {1, 2}], [4, {1, 2, 3}], [6, {1, 2, 3}], [6, {1, 2, 3, 4}], [12, {1, 2, 3, 5}]]

{1, 2, [4, {1, 2}], [4, {1, 2, 3}], [6, {1, 2, 3}], [6, {1, 2, 3, 4}], [12, {1, 2, 3, 5}]}

[infinity, x]

{}

infinity

infinity

table( [ ] )

1

infinity

{}

normal

table( [ ] )

[RootOf(_Z^2+1)/x+1/2, x = x]

{2}

1

{x}

table( [ ] )

0

table( [ ] )

{1}

table( [ ] )

1/2

x

0

1

infinity, 0

{}

normal

[1/2, -1/2, x = x]

0

table( [ ] )

{-1/2, 1/2}

{0}

1

-1/2

1

{1, x}

{-1/2}

{-1}

table( [ ] )

{-2}

table( [ ] )

{-3}

{2}

0

0

{2}

[infinity, 0]

0

[[2, [[[infinity, 1, 1], [0, 1, -1], 0]]]]

2

[]

[2, [[[infinity, 1, 1], [0, 1, -1], 0]]]

[[infinity, 1, 1], [0, 1, -1], 0]

[0, -1/x]

-(1/2)/x

[1, 0, 1]

[0, 4, 0, 1]

[1]

   A Liouvillian solution exists

[4*x^2-1, 4*x, 4*x^2], x, {}, [1], -(1/2)/x, [1, 0, 1], 2, {}, {}, {}, {}, {1, x}

   Group is reducible or imprimitive

[sin(x)/x^(1/2), cos(x)/x^(1/2)]

[sin(x)/x^(1/2), cos(x)/x^(1/2)]

<-- exit \`DEtools/kovacicsols/order2\` (now in \`DEtools/kovacicsols/order2\`) = [sin(x)/x^(1/2), cos(x)/x^(1/2)]}

<-- exit \`DEtools/kovacicsols/order2\` (now in \`DEtools/kovacicsols\`) = [sin(x)/x^(1/2), cos(x)/x^(1/2)]}

<- Kovacic's algorithm successful

y(x) = _C1*sin(x)/x^(1/2)+_C2*cos(x)/x^(1/2)

 

Download dsolve_traced.mw

I suspect that it happens because of the following.

restart;

trace(`simpl/simpl/Re/is`):
#stopat(`simpl/simpl/Re/is`):

Re(x) assuming {};

{--> enter \`simpl/simpl/Re/is\`, args = \`x~\`
                x2 := x

<-- exit \`simpl/simpl/Re/is\` (now in \`simpl/simpl/Re\`) = 0}
                   0

Note that when you enter
    ... assuming {}
the assuming facility temporarily places that property on dependent names. So your call amounts to supplying the property x::{} . So the result is the same as this,

Re(x) assuming x::{};
{--> enter \`simpl/simpl/Re/is\`, args = \`x~\`
                x2 := x

<-- exit \`simpl/simpl/Re/is\` (now in \`simpl/simpl/Re\`) = 0}
                   0

If you follow that procedure in the debugger you can see these checks (which I compare against the case of the unassumed name):

restart;
is(x2,'real');
                 false

restart;
is(x2,'real') assuming x2::{};
                  true

# the previous true induces this test
(not type(x2,'complexcons')
     or `property/ProbablyNonZero`(x2) <> true)
and traperror(is(x2,0)) assuming x2::{};
                  true

# perhaps of particular interest
is(x2,0) assuming x2::{};
                  true

is(x2,0);
                false

You can run it through and see what happens (ie. at lines 4 and 5, say),

showstat(`simpl/simpl/Re/is`);

`simpl/simpl/Re/is` := proc(x)
local x2;
   1   if hastype(x,'float') and not type(x,'complex(numeric)') and type(
         x,'complexcons') then
   2       x2 := `simpl/simpl/ReIm/is/float`(x)
       else
   3       x2 := x
       end if;
   4   if (not type(x2,'complexcons') or has(x2,'RootOf') or 
         `property/ProbablyNonZero`(('Im')(x2)) <> true) and traperror
         (is(x2,'real')) then
   5       if (not type(x2,'complexcons') or has(x2,'RootOf') or 
             `property/ProbablyNonZero`(x2) <> true) and traperror(is(
             x2,0)) then
   6           return 0
           else
   7           return x
           end if
       elif (not type(x2,'complexcons') or has(x2,'RootOf') or 
         `property/ProbablyNonZero`(('Re')(x2)) <> true) and 
         traperror(is(-I*x2,'real')) then
   8       return 0
       end if;
   9   return ('Re')(x)
end proc

Invoking your procedure with the command   wrong();   will only display the one plot that happens to be the return value, since the returned value will be displayed as output.

You could also do this,
   proc() print(plot(x)); printf("z"); print(plot(x)); end proc;
which displays both plots.  (It also returns NULL, since that is the result of the last statement which is a print statement.)

And you can vary it. You could have it forcibly print the first plot and return (and thus incidentally display) the second plot.

Here are some comments, and ways. (I'm sure there are others... possibly simpler).

I chose to use a Vector rather than a vector, and to map the `?[]` indexing command.

vector_of_lists_example_ac.mw

Enter a period (lower dot) for Matrix-Vector and Matrix-Matrix multplication, not a star (which displays as a center-dot in 2D Input).

with(LinearAlgebra):

(Matrix(7, 7, {(1, 1) = -mu, (1, 2) = a[1]*beta[x]*k[1]*PI[h]/mu, (1, 3) = omega[o]*beta[x]*a[1]*PI[h]/mu, (1, 4) = sigma, (1, 5) = -beta[2]*PI[h]/mu, (1, 6) = 0, (1, 7) = a[1]*beta[x]*k[2]*PI[h]/mu, (2, 1) = 0, (2, 2) = -(mu+r[o]+alpha[o])*a[1]*beta[x]*k[1]*PI[h]/mu, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (2, 7) = a[1]*beta[x]*k[2]*PI[h]/mu, (3, 1) = 0, (3, 2) = 0, (3, 3) = -(mu+r[1]+alpha[1]+k[o])*a[1]*beta[x]*k[1]*PI[h]/mu, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = omega[o]*beta[x]*a[1]*k[2]*PI[h]/mu, (4, 1) = 0, (4, 2) = r[o], (4, 3) = r[1], (4, 4) = -mu-sigma, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (5, 1) = 0, (5, 2) = alpha[o], (5, 3) = alpha[1], (5, 4) = 0, (5, 5) = -e[o], (5, 6) = 0, (5, 7) = 0, (6, 1) = 0, (6, 2) = -a[2]*beta[y]*k[1]*PI[m]/mu[b], (6, 3) = -a[2]*beta[y]*k[2]*PI[m]/mu[b], (6, 4) = 0, (6, 5) = -a[2]*beta[y]*PI[m]/mu[b], (6, 6) = -mu[b], (6, 7) = 0, (7, 1) = 0, (7, 2) = a[2]*beta[y]*k[2]*PI[m]/mu[b], (7, 3) = a[2]*beta[y]*k[2]*PI[m]/mu[b], (7, 4) = m[7, 4], (7, 5) = a[2]*beta[y]*PI[m]/mu[b], (7, 6) = 0, (7, 7) = -mu[b]})).(Matrix(7, 1, {(1, 1) = omega[1], (2, 1) = omega[2], (3, 1) = omega[3], (4, 1) = omega[4], (5, 1) = omega[5], (6, 1) = omega[6], (7, 1) = omega[7]})) = (Matrix(7, 1, {(1, 1) = 0, (2, 1) = 0, (3, 1) = 0, (4, 1) = 0, (5, 1) = 0, (6, 1) = 0, (7, 1) = 0}))

(Matrix(7, 1, {(1, 1) = -mu*omega[1]+a[1]*beta[x]*k[1]*PI[h]*omega[2]/mu+omega[o]*beta[x]*a[1]*PI[h]*omega[3]/mu+sigma*omega[4]-beta[2]*PI[h]*omega[5]/mu+a[1]*beta[x]*k[2]*PI[h]*omega[7]/mu, (2, 1) = -(mu+r[o]+alpha[o])*a[1]*beta[x]*k[1]*PI[h]*omega[2]/mu+a[1]*beta[x]*k[2]*PI[h]*omega[7]/mu, (3, 1) = -(mu+r[1]+alpha[1]+k[o])*a[1]*beta[x]*k[1]*PI[h]*omega[3]/mu+omega[o]*beta[x]*a[1]*k[2]*PI[h]*omega[7]/mu, (4, 1) = r[o]*omega[2]+r[1]*omega[3]+(-mu-sigma)*omega[4], (5, 1) = alpha[1]*omega[3]+alpha[o]*omega[2]-e[o]*omega[5], (6, 1) = -a[2]*beta[y]*k[1]*PI[m]*omega[2]/mu[b]-a[2]*beta[y]*k[2]*PI[m]*omega[3]/mu[b]-a[2]*beta[y]*PI[m]*omega[5]/mu[b]-mu[b]*omega[6], (7, 1) = a[2]*beta[y]*k[2]*PI[m]*omega[2]/mu[b]+a[2]*beta[y]*k[2]*PI[m]*omega[3]/mu[b]+m[7, 4]*omega[4]+a[2]*beta[y]*PI[m]*omega[5]/mu[b]-mu[b]*omega[7]})) = (Matrix(7, 1, {(1, 1) = 0, (2, 1) = 0, (3, 1) = 0, (4, 1) = 0, (5, 1) = 0, (6, 1) = 0, (7, 1) = 0}))

NULL


 

Download trc_ac.mw

I found that using solve might depend on whether the changed equation was expanded, and upon the working precision (Digits).

You can also adjust Digits to see effect on residual (forward) error.

There are other ways to compute this (eg. using procedures which each raise working precision ie. Digits locally while keep fsolve's notion of accuracy (dependent on inbound Digits) coarser. I haven't done that. Similar effects can be had by treating it as an optimization problem (eg. by using DirectSearch, or setting it up that way).

restart;

with(LinearAlgebra):

TrainLoad := -10*10^6*(cos(convert(40*degrees, radians))+I*sin(convert(40*degrees, radians)));

-10000000*cos((2/9)*Pi)-(10000000*I)*sin((2/9)*Pi)

 

evalf(TrainLoad, 7);

-7660444.-6427876.*I

f1n2 := (0.03 + I*0.1515)*Ix[c1] - (0.03 + I*0.1515)*Ix[c2] + 2 * V[at1] - 55*10^3 = 0;

(0.3e-1+.1515*I)*Ix[c1]+(-0.3e-1-.1515*I)*Ix[c2]+2*V[at1]-55000 = 0

f3n4 := (1.6 + I*6.24)*Ix[c1] + (1.12 + I*2.64)*Ix[c2] + V[t] - V[at1] = 0;

(1.6+6.24*I)*Ix[c1]+(1.12+2.64*I)*Ix[c2]+V[t]-V[at1] = 0

f5n6 := (1.36 + I*4.44)*Ix[c2] + V[at2] - V[t] = 0;

(1.36+4.44*I)*Ix[c2]+V[at2]-V[t] = 0

f7n8 := (-1.12 - I*2.64)*Ix[c1] + (-3.92 - I*12.00)*Ix[c2] + V[at2] - V[at1] = 0;

(-1.12-2.64*I)*Ix[c1]+(-3.92-12.00*I)*Ix[c2]+V[at2]-V[at1] = 0

f9n10 := V[t] * (Ix[c1] - Ix[c2]) + TrainLoad = 0;

V[t]*(Ix[c1]-Ix[c2])-10000000*cos((2/9)*Pi)-(10000000*I)*sin((2/9)*Pi) = 0

polynomials := {f1n2, f3n4, f5n6, f7n8, f9n10};

{(1.36+4.44*I)*Ix[c2]+V[at2]-V[t] = 0, V[t]*(Ix[c1]-Ix[c2])-10000000*cos((2/9)*Pi)-(10000000*I)*sin((2/9)*Pi) = 0, (-1.12-2.64*I)*Ix[c1]+(-3.92-12.00*I)*Ix[c2]+V[at2]-V[at1] = 0, (0.3e-1+.1515*I)*Ix[c1]+(-0.3e-1-.1515*I)*Ix[c2]+2*V[at1]-55000 = 0, (1.6+6.24*I)*Ix[c1]+(1.12+2.64*I)*Ix[c2]+V[t]-V[at1] = 0}

variables := {Ix[c1], Ix[c2], V[at1], V[at2], V[t]};

{Ix[c1], Ix[c2], V[at1], V[at2], V[t]}

fsolve(polynomials, variables, complex);

{Ix[c1] = 955.2297105-5281.491898*I, Ix[c2] = -505.0156845+2424.830843*I, V[at1] = 26894.34238+4.981252431*I, V[at2] = 10829.70666+56.66545127*I, V[t] = -623.3636109+1112.165758*I}

f9n10 := V[t] * conjugate(Ix[c1] - Ix[c2]) + TrainLoad = 0;

V[t]*conjugate(Ix[c1]-Ix[c2])-10000000*cos((2/9)*Pi)-(10000000*I)*sin((2/9)*Pi) = 0

Digits:=20:

sol1:=solve({f1n2, f3n4, f5n6, f7n8, f9n10}, variables);

Warning, solutions may have been lost

sol1:=solve({f1n2, f3n4, f5n6, f7n8, expand(f9n10)}, variables);

{Ix[c1] = 1236.8820951263380105-4880.5299795094608907*I, Ix[c2] = -630.06302807299903300+2236.4267367463552561*I, V[at1] = 26932.886351895631871-34.66674233851253885*I, V[at2] = 11895.825533339687912-1029.4551170858521807*I, V[t] = 1109.205104006591891-785.3945997549247388*I}

K:=map(var->var=cat(re,var)+I*cat(im,var),variables):
evalc(eval(map(eq->[Re(eq),Im(eq)][],{f1n2, f3n4, f5n6, f7n8, f9n10}),K)):
sol2:=eval(K,simplify(fnormal(fsolve(%,complex)),zero));

{Ix[c1] = 1236.8820951263380105-4880.5299795094608906*I, Ix[c2] = -630.06302807299903300+2236.4267367463552561*I, V[at1] = 26932.886351895631871-34.666742338512538846*I, V[at2] = 11895.825533339687912-1029.4551170858521808*I, V[t] = 1109.2051040065918903-785.39459975492473915*I}

evalf[100](eval((lhs-rhs)~({f1n2, f3n4, f5n6, f7n8, f9n10}), sol1)):
max(abs~(evalf[5](%)));

0.56912801178645213474e-11

evalf[100](eval((lhs-rhs)~({f1n2, f3n4, f5n6, f7n8, f9n10}), sol2)):
max(abs~(evalf[5](%)));

0.16209122737779487861e-12

 

Download no_conjugate_ac.mw

You can also use the command-completion facility, or the palettes, or the Ctl-Shift-G keyboard shortcut to easily construct compound names in 2D Input mode.

This also allows the names to appear in their nice typeset form as 2D Input, and not just as 2D Output.

See Help pages for command-completion and the entry palettes.

restart;

 

For the following 2D Input I typed in the keystrokes D e l t a and

then I pressed the Escape key to activate command-completion .

 

That changed it to the typeset Greek Delta, and then I finished

entering rest of the name.

 

`&Delta;x`

`&Delta;x`

 

You can see the underlying structure of the name by using the

lprint command.

 

lprint(%);

`&Delta;x`

 

Below, I entered the first portion using the Greek entry palette

(from the GUI's left-panel).

 

`&Delta;x`

`&Delta;x`

lprint(%);

`&Delta;x`

 

Download cmd_completion.mw

As Carl mentioned, you can also enter the compound names using its 1D equivalent notation (which -- as you can see above -- involves a name containing HTML-style entity, here). But for some kinds of typeset name you may not find that easy to guess, eg. names with a caret over them, etc.

In 2D Input mode, you can also use the keyboard acceleration Ctl-Shift-G to temporarily place the GUI into Greek-entry-mode. For your example, the same result for typeset 2D Input can also be had using the keystrokes,
   Ctl-Shift-G D x
There, the "D" makes it enter a Greek uppercase Delta. See the Help page.

In a previous Question you were shown in an Answer that the simplify command accepts the assume=[...] option.

That does not mean that option assume=[...] works for all commands. It does not.

The odetest command does not throw an error if you pass it some extra arguments.

If you want to construct your call (to utilize the assuming facility) from a precomputed list or collection of inputs then you can call using its prefix form:

sol:=y(x)^(1/2)=-1+2*exp(x):
ode:=diff(y(x),x)-2*y(x)=2*y(x)^(1/2):
`assuming`([odetest(sol,ode)], [x::real, x>0]);
                     0

Why are you trying to pass the plot function call around as if it were an object!? Do you have some special purpose where there really isn't a better programming methodology? (I can think of very few, and since you are new to Maple I suspect that your actual goal could be better done in another way altogether).

Having said that,

  eval('plot'(x, x = 0 .. a), a = 3);
  value(eval(%plot(x, x = 0 .. a), a = 3));
  ((e::uneval)->eval(e,a=3))( plot(x, x = 0 .. a) );

On this forum it is almost always better to inform us in detail of your actual motivating goal, rather than just inquiring about some barrier found during an attempt you devised for getting there.
 

You can customize and adjust this.

The idea is that you can write a convenient procedure to generate such plots in a reusable manner. And you could also augment the results futher.

restart;

F:=proc(ee,var,a,b,n::posint,
        {lineoptions::list:=[],
         ptoptions::list:=[]})
  local i,pts,xpts;
  uses plots;
  xpts:=[seq(a+(i-1)*(b-a)/(n-1),i=1..n)];
  try
    pts:=map(p->[p,eval(ee,var=p)],xpts);
  catch:
    pts:=map(p->[p,limit(ee,var=p)],xpts);
  end try;
  display(plot(ee,var=a..b,_rest),
          seq(plottools:-line([p[1],0],p,lineoptions[]),p=pts),
          pointplot(pts,ptoptions[]),
          ':-axis[1]'=[':-tickmarks'
                       =(u->u=Typesetting:-Typeset(u))~(xpts)]);
end proc:

F(sin(x),x,-2*Pi,2*Pi,7);

F(cos(x),x,-2*Pi,2*Pi,9,color=green,
  lineoptions=[color=blue,linestyle=dash],
  ptoptions=[color=red,symbol=solidcircle,symbolsize=10],
  size=[400,300]);

F(sin(x)/x,x,-2*Pi,2*Pi,9,color=blue,thickness=3,
  lineoptions=[color=blue,thickness=3],
  ptoptions=[color=blue,symbol=solidcircle,symbolsize=20],
  size=[500,200], labels=["",""]);

plots:-display(
  F(sin(x)/x,x,-2*Pi,2*Pi,13,color=blue,thickness=3,
    lineoptions=[color=blue,thickness=3],
    ptoptions=[color=blue,symbol=solidcircle,symbolsize=20],
    size=[600,300], labels=["",""], axes=box),
  plot(0,-2*Pi..2*Pi,color=black),
  plottools:-line([0,1],[0,-0.3],color=black),
  axis[2]=[tickmarks=[0,0.5,1.0]]
);

 

Download plotting_fun.mw

restart;

with(Units:-Standard):

a := 153350 * Unit(J);

153350*Units:-Unit(J)

#
# This doesn't work. (It's a mistaken attempt.)
#
convert(a,'units',Unit('kN')*Unit('m'));

(3067/20)*Units:-Unit(1000*J)

#
# This works
#
convert(a,'units',Unit('kN'*'m'));

(3067/20)*Units:-Unit(kN*m)

#
# The first attempt is a mistake, because the
# 3rd argument is itself combined, unitwise.
#
Unit('kN')*Unit('m');

1000*Units:-Unit(J)

#
# Alternatively, you can just use the synbols.
#
convert(a,'units','kN'*'m');

(3067/20)*Units:-Unit(kN*m)

It can also be done in 2D Input mode.

convert(a, 'units', 'kN'*'m')

(3067/20)*Units:-Unit(kN*m)

It can also be done sing the Units palette, by choosing
the blue "units" item. Hence the units appear in upright
Roman in Maple 2020 (or braced, when inlined on Mapleprimes).

convert(a, 'units', Unit('kN'*'m'))

(3067/20)*Units:-Unit(kN*m)

 

Download unitconvert.mw

Since NULL is ephemeral when used in ways that you intend (eg. being passed around, or returned, in an expression sequence) then it isn't a very good choice. Don't choose something that can be problematic for one of your basic usage patterns.

Using {} is not terrible, but it will fall down for a computation where {} is also a possible natural result and where you need to distinguish between the two situations. In that case you'd have to utilize something else instead, which means that your methodologies could be inconsistent. Note that you didn't cite only computations where the result denoted possible solutions -- so {} might make sense to denote the empty set, but is not great to denote general cases where something "could not be computed".

How about using FAIL instead? Note that your two cited example circumstances were, "...it could be something that could not be computed inside foo, or some solution that could not be found". The common phrasing there is "could not be".

For the cases of computing possible solutions, would you want to be able to distinguish between the case of provably no solutions and the case where it happens to not succeed in finding any solution?

[edit] I don't think that solutions involving returning uneval-quoted NULL is great. You'd need to avoid using whatever it got assigned to in an analogous situation. Why introduce any need for extra care, if you don't have to?

Here is one way to deal with that, which consists of the step,

   subsindets(eq[4], float,u->`if`(frac(u)=0.0,trunc(u), u));

(There are more complicated ways, some of which deal with other things that you might encounter and want, such as stripping off trailing zeroes. Let us know...)

restart

with(LinearAlgebra):

par := {a = 2.5, alpha = 2, k_r = .5, k_rc = .2, k_rq = .2, kappa = 0.1e-2, mu_k = .35, mu_s = .44, omega = .766620580157922, sigma_0 = 110, sigma_1 = 1.37, sigma_2 = 0.823e-1, x_s3 = -1.04003626422324936017819852700633621040584050594846801927800, zeta = 0.904504977553123318334601762827181333680096702957228781770315e-1}:

f := proc (v_r) options operator, arrow; mu_k+(mu_s-mu_k)*exp(-a*v_r) end proc:

``

g_exp1 := taylor(1/f(v_rv+x21), x21 = 0, 4):

for k to 4 do g_coeff[k] := taylor(subs(v_r = v_rv, coeff(g_exp1, x21, k-1)), epsilon = 0, 3) end do:

g0 := eval(subs(par, coeff(g_coeff[1], epsilon, 0))):

g1 := eval(subs(par, coeff(g_coeff[2], epsilon, 0))):

g2 := eval(subs(par, coeff(g_coeff[3], epsilon, 0))):

g3 := eval(subs(par, coeff(g_coeff[4], epsilon, 0))):

``

NULL

eq[1] := subs(par, diff(x[1](t), t)-x[2](t)):

eq[2] := subs(par, diff(x[2](t), t)+2*zeta*x[2](t)+x[1](t)+k_r*(x[1](t)-x[3](t))+2*kappa*(x[2](t)-x[4](t))+k_rq*(x[1](t)-x[3](t))^2+k_rc*(x[1](t)-x[3](t))^3)

diff(x[2](t), t)+.1829009955*x[2](t)+1.5*x[1](t)-.5*x[3](t)-0.2e-2*x[4](t)+.2*(x[1](t)-x[3](t))^2+.2*(x[1](t)-x[3](t))^3

eq[3] := subs(par, diff(x[3](t), t)-x[4](t))

diff(x[3](t), t)-x[4](t)

eq[4] := subs(par, diff(x[4](t), t)+2*kappa*alpha*(x[4](t)-x[2](t))+k_r*alpha*(x[3](t)-x[1](t))-k_rq*alpha*(x[3](t)-x[1](t))^2+k_rc*alpha*(x[3](t)-x[1](t))^3+alpha*(sigma_0*x[5](t)+sigma_1*v_r*(1-sigma_0*x[5](t)*(g_0+g_1*x[4](t)+g_2*x[4](t)^2+g_3*x[4](t)^3))+sigma_2*v_r)):

diff(x[4](t), t)+0.4e-2*x[4](t)-0.4e-2*x[2](t)+x[3](t)-x[1](t)-.4*(x[3](t)-x[1](t))^2+.4*(x[3](t)-x[1](t))^3+220*x[5](t)+2.74*v_r*(1-110*x[5](t)*(g_0+g_1*x[4](t)+g_2*x[4](t)^2+g_3*x[4](t)^3))+.1646*v_r

eq[5] := subs(par, diff(x[5](t), t)-v_r*(1-sigma_0*x[5](t)*(g_0+g_1*x[4](t)+g_2*x[4](t)^2+g_3*x[4](t)^3)))

diff(x[5](t), t)-v_r*(1-110*x[5](t)*(g_0+g_1*x[4](t)+g_2*x[4](t)^2+g_3*x[4](t)^3))

NULL

Download question_sunit_ac1.mw

Here is one way.

Let me know if you really prefer that multiple calls foo(ode,y,x) would produce the exact same result (in which case it could adjusted to be store and manage that, I think.)

restart;

foo := proc(ode::`=`,y,x::symbol)
  local t:=`tools/genglobal`(x);
  if t=x then t:=`tools/genglobal`(x); end if;
  return PDEtools:-dchange({x=t^2},ode,t)  
end proc:

ode:=diff(y(x),x)=sin(x);
foo(ode,y,x);
foo(ode,y,x);

diff(y(x), x) = sin(x)

(1/2)*(diff(y(x0), x0))/x0 = sin(x0^2)

(1/2)*(diff(y(x1), x1))/x1 = sin(x1^2)

ode:=diff(y(t),t)=sin(t);
foo(ode,y,t)

diff(y(t), t) = sin(t)

(1/2)*(diff(y(t0), t0))/t0 = sin(t0^2)

 

Download genglobal.mw

First 113 114 115 116 117 118 119 Last Page 115 of 336