Joe Riel

9530 Reputation

23 Badges

20 years, 24 days

MaplePrimes Activity


These are answers submitted by Joe Riel

Here are a few suggestions. First, assign N directly and then use it to compute dt. Second, use the evalf to force Maple to convert 2*Pi/N to a floating point. Change vector to Vector. The function assignment should be
Lambda := (t,x,y) -> sqrt((cos(t)-x)^2+(sin(t)-y)^2);
With those changes everything seemed to work.
One way is to use the save and read commands:
y := x^2:
save y, "somefile.mpl":
restart;
read "somefile.mpl":
y;
              x^2
Use fsolve. For example,
deq := diff(y(t),t)=-y(t):
dsol1 := dsolve({deq, y(0)=1}, numeric, output=listprocedure):
dsol2 := dsolve({deq, y(0)=2}, numeric, output=listprocedure):
delta := eval(y(t), dsol2) - eval(y(t), dsol1): # I swapped order so delta is positive
fsolve(delta=0.5);
                  0.6931469270
delta(%); # check result
                   0.5000000000
Note that you can pass an interval to fsolve to restrict its search to a specified real interval. See ?fsolve,details.
Here's one approach:
deq := diff(y(t),t)=-y(t):
dsol1 := dsolve({deq, y(0)=1}, numeric, output=listprocedure):
dsol2 := dsolve({deq, y(0)=2}, numeric, output=listprocedure):
delta := eval(y(t), dsol1) - eval(y(t), dsol2):
plot(delta, 0..10);
First, traperror is obsolete, the proper way to catch errors is with try/catch. A simple way to "trap" warnings is to replace the WARNINGS procedure. Here is code I use to capture all warnings when evaluating a particular expression. The try statement is necessary to handle errors that occur in the evaluation (this is used during testing).
GetWarning := proc(expr::uneval)
global WARNING;
local warnings,oldWARNING,cnt,i;
description "Capture arguments passed to calls to WARNING."
    ,"Return a list of lists, each sublist corresponds to the arguments "
    "to one call to WARNING.";
    unprotect(WARNING);
    oldWARNING := eval(WARNING);
    cnt := 0;
    warnings := table();
    try
        WARNING := proc()
            cnt := cnt+1;
            warnings[cnt] := [args];
            NULL;
        end proc;
        eval(expr);
    finally
        WARNING := eval(oldWARNING);
        protect('WARNING');
    end try;
    [seq(warnings[i], i=1..cnt)];
end proc:

f := proc()
    WARNING("danger!");
    WARNING("alert!");
end proc:

GetWarning(f());
               [["danger!"], ["alert!"]]
What does the parenthetical phrase, "(xmaple -classic)," mean? Do you launch it from a shell prompt with that command? On my debian box, that launches the standard interface, not the classic (the -classic option does nothing). To start the classic interface I do
maple -cw
To clarify Jacques' warning about positioning, consider the ODE
sys := {diff(x(t),t)=sin(t),diff(y(t),t)=cos(t),x(0)=1,y(0)=1}:
and its numeric solution
dsol := dsolve(sys, numeric):
dsol(1); 
          [t = 1., x(t) = 1.45969768100667418, y(t) = 1.84147131233092276]
You do not want to assume that x(t) will always be the second element in the list (for a particular solution, the order won't change, but a subsequent call to dsolve could conceivably change the order of x(t) and y(t), depending on how the ordering was created). To create a function of x(t) and y(t) you might do
S := tt -> eval(3*x(t) + 4*y(t), dsol(tt)):
The usual way to plot an output of dsolve/numeric is to use plots[odeplot]:
plots[odeplot](dsol, [t, 3*x(t)+4*y(t)], 0..1)
If I understand correctly, there isn't a way to do what you want. That is, given an existing module there is no straightforward way to add an export to that module. Even modifying an export of a module isn't generally possible (there are some hackish ways, but they have limits). This is a drawback of the modular approach compared to the old table-based package approach used in Maple R5—while the modular approach is more robust and better from the programmar's (creator's) viewpoint, it is less malleable for the user. It is, of course, possible to add a new module to an existing archive, though usually it is better to just create a new archive and insert it in the libname path. Use LibraryTools[Save]. To make a module ``withable,'' give it the package option. For example:
mymodule := module()
export A,B:
option package;
  A := proc() ... end proc:
  B := proc() ... end proc:
end module:
LibraryTools:-Save(mymodule,"mylib.mla"):
The ability to rebind operators, such as +, was introduced in Maple 6; however, the overload option (and procedure) was introduced in Maple 9.5.
See the help page for NextAfter.
> Digits := 20:             
> epsilon := NextAfter(0,1);
                                                                   -20
                                                  epsilon := 0.1 10

> 1+epsilon;                
                                                 1.0000000000000000000
That isn't quite what you want, but...
I'm not sure I understand, but the following (rather ugly) code does what I think you want
restart; with(stats); Digits := 14;

N := 4:
L := 0.002:
beta:=0.85:
cost := Vector([200,100,700,1000]):
lambda := Vector([0.2,0.6,0.3,0.7]):
Lambda := add(lambda[i], i=1..N):

DFR := s -> ( add(lambda[i]*add((L*lambda[i])^k*exp(-L*lambda[i])/k!
                                , k = 0 .. s[i]-1)
                  , i = 1 .. N)
              / Lambda );

IN := Matrix(N,N,shape=identity):
S := Vector[row](4):
maxval := DFR(S);

while maxval < beta do
    prevval := maxval;
    maxrelval := 0;
    for k to N do
        Sk := S + IN[k,1..N];
        val := DFR(Sk);
        relval := (val - prevval)/cost[k];
        if relval > maxrelval then
            maxrelval := relval;
            maxval := val;
            Snxt := Sk;
        end if;
    end do;
    if maxrelval = 0 then
        print("cannot increase");
        break;
    end if;
    S := Snxt;
end do:

S, maxval;
I think the following does what you want.
restart; with(stats); Digits := 14;

N := 4:
L := 0.002:
beta:=0.85:
lambda := Vector([0.2,0.6,0.3,0.7]):
Lambda := add(lambda[i], i=1..N):

DFR := s -> ( add(lambda[i]*add((L*lambda[i])^k*exp(-L*lambda[i])/k!
                                , k = 0 .. s[i]-1)
                  , i = 1 .. N)
              / Lambda );

IN := Matrix(N,N,shape=identity):
S := Vector[row](4):
maxval := DFR(S);

while maxval < beta do
    found := false; # flag to avoid infinite loop if no Sk is better than S.
    for k to N do
        Sk := S + IN[k,1..N];
        val := DFR(Sk);
        if val > maxval then
            maxval := val;
            Snxt := Sk;
            found := true;
        end if;
    end do;
    if not found then break; end if; # may want to add warning
    S := Snxt;
end do;

maxval;
I don't understand what you are attempting, however, there seem to be two clear problems. First, you probably want to use add instead of sum in the computation of DFR (you used one of each). Second, your do loop is confusing. It's index is i, but you also use i as the index in the adds. These are different i's and I'm not sure which you intend. Change one of them to a j to make it clear what you are doing.
A good method for creating types is to use the TypeTools module, in particular, its AddType export. You'd do something like:
Automatism := module()
export zeros;
local ModuleLoad,ModuleUnload;

  ModuleLoad := proc()
    TypeTools:-AddType( 'TransferFunction', handler );
  end proc;

  ModuleUnload := proc()
    TypeTools:-RemoveType( 'TransferFunction' );
  end proc;

  zeros := proc( f:: TransferFunction) ... end proc;

  ModuleLoad();
end module:
You have to write the 'handler'; see the AddType help pages for details (it wasn't clear to me what your data structure is, otherwise I'd have suggested something). ,
You'd use seq something like this:
display([ seq( plot([X1[k],Y1[k]],[X2[k],Y2[k]]), k=1..20) ], 'insequence');
First 103 104 105 106 107 108 109 Last Page 105 of 114