Joe Riel

9660 Reputation

23 Badges

20 years, 4 days

MaplePrimes Activity


These are answers submitted by Joe Riel

This is known as the "money-changing problem".  See Wilf's "GeneratingFunctionology" for a discussion of solving it with generating functions (per Robert's approach). You can download an older edition of Wilf's book from his site, but I recommend purchasing the book, it is both informative and enjoyable.

By default, prime notation represents differention with respect to x, not t. You want to use the dot notation.  To do that, select the symbol, type Ctrl-Shift-" (the last charater is a double-quote) then type the appropriate number of dots (periods).

How do you plan to use this model?  That is, is it part of a larger simulation? 

One way to model simple hysteresis is with events in dsolve, see ?dsolve/events.  For example

evts := [ [ i(t) - ihys(t), ihys(t)=-ihys(t)]  ]:
dsys := { diff(i(t),t)=sin(t), i(0)=-1, ihys(0)=1/2 }:
integ := dsolve(dsys, numeric
                , 'events' = evts
                , 'discrete_variables' = [ ihys(t) ]
               ):
plots:-odeplot(integ, [[t,i(t)],[t,ihys(t)]]);

Modeling the continuous hysteresis of a B-H loop in ferromagnetic material is more challenging. A typical approach is to use the Jiles-Atherton equations

Specify the colors as a list:   color=[blue,green,red].  

Attach a "World Applied Force" block to the end and set its resolved frame to "outboard".  That sets the frame of the resolved force to that of the frame to which it is attached.  Then drive the z-input of the force block with the desired orthogonal force.

An amusing way, which mostly works, is

eval(eval({op(Ex)},G=(x->exp(x))),exp=1);

But y=0 does not solve the equation. 

 

Peter's method is nice.  Less nice, here, but useful in other applications, is to use the two-argument form of ?op.  The first argument specifies the operand, it may be an integer, a range, or a list of integers, with the final element allowed to be a range.  Here you would do

op([1,1,1,2],minimum[2])

There is an issue with the pseudo-code: the internal summation uses the index i, which is also the index of the outer loop.  That cannot work properly, however, it is not entirely clear what is intended.

FPackage :=  module()
option package;
export `*`;
    `*` := overload([NULL
                     , proc(a::algebraic, b::specfunc(anything, f))
                       option overload;
                           map2(procname, a, b);
                       end proc      
                     , proc(b::specfunc(anything, f), a::algebraic)
                       option overload;
                           map2(procname, a, b);
                       end proc
                    ]
                   );
end module:

with(FPackage);
a*f(b,c,d,f(e,g));
                               f(a*b,a*c,a*d,f(a*e,a*g))

One way is to delay the evaluation of the function call norma by quoting it with single forward quotes:

int( 'norma'(...), ..., numeric )

 

A slightly different approach is to use ?solve/identity

ex := R*d^2 / ( R*d^2 + s*L + s^2*R*L*C ):

solve(identity(ex = 1/(1 + 2*zeta*(s/omega)  + (s/omega0)^2), s), {zeta,omega0});
               d             omega L                  d             omega L
  {omega0 = --------, zeta = -------}, {omega0 = - --------, zeta = -------}
                 1/2            2                       1/2            2
            (L C)            2 d  R                (L C)            2 d  R

y := 23.45*sin((360*(x+284))/284):
seq('y'(x) = evalf(y), x = 1..365);

Create a system with both systems.  For example

with(DynamicSystems):
sys1 := TransferFunction(1/(s+2)):
sys2 := TransferFunction(s/(s+3)):
sys12 := TransferFunction(<sys1:-tf,sys2:-tf>):
BodePlot(sys12);

 



Use a ?piecewise expression to represent y in terms of x.  Thus

y := piecewise(x<=4, 9/10, 1-1/10*exp(-x/10)):

You can integrate and plot this expression as usual.

First 49 50 51 52 53 54 55 Last Page 51 of 114