Carl Love

Carl Love

28055 Reputation

25 Badges

12 years, 357 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Here's another recursive procedure that solves the problem.

bit:= proc(n::nonnegint, k::nonnegint)
local
   Put:= proc(T,n,v) T[n]:= v; T end proc,
   recur:= proc(n,k)
   local j; 
      if n=k then {table([seq(j=1, j= 1..n)])}
      elif k=0 then {table([seq(j=0, j= 1..n)])}
      else map(Put, thisproc(n-1, k), n, 0) union map(Put, thisproc(n-1, k-1), n, 1)
      end if
   end proc
;
   convert~(recur(n,k), list)
end proc:       

 

The time and memory values on the status line are only updated after a garbage collection. If the computation is at a stage where no garbage is being generated, then there'll be no garbage collection and no status updates.

If you're not familiar with the computer-science lingo, garbage is memory that is no longer accessible. For example, if you did

A:= [1,2];  A:= [3,4];

then the memory holding [1,2] would become garbage. Garbage collection is the process of collecting these scraps of unused memory and rearranging them into usable chunks of available memory. Garbage collection only occurs after a certain amount of garbage has piled up, typically many megabytes.

When you replace a list element with NULL, it changes the length of the list:  it is reduced by one. This behavior is different for Vectors: NULL can be a Vector element. 

I think that you can achieve what you want without a lot of intricate work with tildes. The procedure Vectorize, below, is essentially a form of unapply that takes a symbolic procedure and returns a procedure that's designed to work efficiently on float[8] rtables but will also work on other rtables or container objects. I believe that this completely automates the operation that you described.

Vectorize:= P-> proc(X) try map[evalhf](P,X) catch: P~(X) end try end proc:
A:= LinearAlgebra:-RandomMatrix((3000,3000), datatype= float[8]):
F:= Vectorize(D(x-> sin(x^2))):
FA:= CodeTools:-Usage(F(A));

memory used=0.54GiB, alloc change=68.67MiB, cpu time=4.33s, real time=3.67s, gc time=1.48s

 

The corresponding input is Eval(H1, t= 0.45). That's with a capital E so that it's inert. To activate the evaluation, you'd use value(%) (if you were actually running Maple, of course).

Use the -q (for quiet) flag on the cmaple command, for example

cmaple -q < myinputfile.mpl

In your module, include a procedure named ModuleLoad (you must use precisely this name). It may be declared local or export, but it needs to be one of these. Inside this procedure, include the statement :-Clr:= B. Right before the end module, include the statement ModuleLoad().

I don't understand why you want Clr to be global rather than a module export.

You need to use eval rather than subs. It is a lot like subs, but much more mathematically sophisticated, and it understands how to handle derivatives, integrals, etc. In constrat, subs is intended primarily for low-level programming. It just blindly substitutes what you tell it to.

Here's how to use eval in your situation:

g:= (x,y)-> diff(u1(x, y), x, x) + diff(u2(x, y), x, y):
f:= (y)-> eval(g(x,Y), [x= 0, Y= y]);

Note the distinction that I have made between lowercase y and uppercase Y.

Here's a completely different way:

g:= D[1,1](u1) + D[1,2](u2):
f:= (y)-> g(0,y):

 

 

 

Here's an elegant and automated solution to your problem. The main idea is to create a new ODE and a new BC by differentiating the expression that you want evaluated.
 

An iterative solution to a parametrized BVP

Carl J Love 2016-Dec-22


restart:


#Your original system:
#(This should be the only thing you that need to modify):

EQs:=
   diff(f(x),x$3) + f(x)*diff(f(x),x$2) - a*diff(f(x),x$1)^2 = 0,
   diff(g(x),x$2) + b*f(x)*diff(g(x),x$1) = 0
;
LB:= 0:  UB:= 5: #upper and lower boundaries
BCs:= f(LB)=0, D(f)(LB)=1, D(f)(UB)=0, g(LB)=0.5, g(UB)=0;
expr:= a*f(x) + b*diff(g(x),x) + c*diff(f(x),x)*g(x);
#expr to be evaluated at:
x__e:= 1.; Params:= [a,b,c]=~ ([seq(0..1, 0.2)] $ 3);

diff(diff(diff(f(x), x), x), x)+f(x)*(diff(diff(f(x), x), x))-a*(diff(f(x), x))^2 = 0, diff(diff(g(x), x), x)+b*f(x)*(diff(g(x), x)) = 0

f(0) = 0, (D(f))(0) = 1, (D(f))(5) = 0, g(0) = .5, g(5) = 0

a*f(x)+b*(diff(g(x), x))+c*(diff(f(x), x))*g(x)

1.

[a = [0, .2, .4, .6, .8, 1.0], b = [0, .2, .4, .6, .8, 1.0], c = [0, .2, .4, .6, .8, 1.0]]


#Implicity define a new function h(x) to be expr, and make it into
#a new ODE:
Expr:= h(x) = expr:
newEQ:= diff(h(x) = Expr, x);

#Create a BC for h(0) by evaluating Expr at x=0 and then evaluating
#wrt BCs:
#(This step can be tricky; may be difficult to fully automate.)
#(Might also want to try to make a new BC at UB if the one at LB
#doesn't work.)
newBC:= eval(convert(eval(Expr, {x= LB}), D), {BCs});

#Abstract the BVP system parametrized by (a,b,c):
sys:= unapply({EQs, newEQ, BCs, newBC}, [a,b,c]):

#And use that to create a parametrized solver that explicitly
#evaluates h(x) at a point.
_H:= (a::realcons, b::realcons, c::realcons, x0::realcons)->
   eval(h('x'), dsolve(sys(a,b,c), numeric)(x0))
:

#Example usage:
h(x__e) = _H(0.5, 0.5, 0.5, x__e);
 

diff(h(x), x) = (diff(h(x), x) = a*(diff(f(x), x))+b*(diff(diff(g(x), x), x))+c*(diff(diff(f(x), x), x))*g(x)+c*(diff(f(x), x))*(diff(g(x), x)))

h(0) = b*(D(g))(0)+.5*c

h(1.) = HFloat(0.31148732145735375)

#Do calculations, putting them into a table:
Calc:= proc()
local abc, ABC, Htable:= table();

   for abc in Iterator:-CartesianProduct(eval([a,b,c], Params)[]) do
      ABC:= convert(abc, list);
      assign(`?[]`(Htable, ABC) = _H(ABC[], x__e))
   end do;
   Htable
end proc:

#Measure resource consumption for the "hard" computation:
Htable:= CodeTools:-Usage(Calc()):

memory used=431.50MiB, alloc change=0 bytes, cpu time=5.44s, real time=5.46s, gc time=218.75ms


#Put the table in an elegantly displayed format:
N:= mul(nops~(eval([a,b,c], Params))): #number of parameter tuples.
sav:= interface(rtablesize= N+2):
DataFrame(
   Matrix((`[]`@op)~([entries(Htable, pairs, indexorder)])),
   columns= [a, b, c, h(x__e)], rows= ['convert(``, `local`)' $ N]
);
interface(rtablesize= sav):

RTABLE(18446744074183065598, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 217, 1 .. 5)

 


 

Download ParamtrizedBVP.mw

A direct translation into Maple of what you have written is

local D:= n->
   Matrix(
      ((n+1)$2),
      (i,k)->
         if i <> k then P[n](t[i])/P[n](t[k])/(t[i] - t[k])
         elif i = 1 then -n*(n+1)/2
         elif i = N+1 then n*(n+1)/2
         else 0
         end if
   )
:

I can't read your denominators on n*(n+1). I have guessed 2.

Your method of presentation shows disrespect for this forum and its members. If you can't come up with a better presentation, you should at least first ask for permission to post something so crude, and beg for the forgiveness of your readers.

When you make an assignment to F[1], it also changes the meaning of F, which becomes a table, not a vector field. If you want a new variable which displays as F with a subscript, use F__1 (that's two underscores). That is a completely new variable, with no connection to F. When any name is followed by [anything], the resulting structure has some relation to the parent name.

Regarding your title: Consider this simple example:

x:= 2:
x:= 1/(x-1);
x:= 1/(x-1);

Surely you can understand that it is correct that the second command works and that the third one returns an error although they are syntacticly the same command.

One way to do what you want (and it's certainly the easiest programmatic way) is make a Matrix whose first row is the headers and whose second row is the matrices. So, some of the elements of this Matrix are themselves matrices.

Display:= Matrix(2,3, [[`Matrix 1`, `Matrix 2`, `Matrix 3`]]):
Display[2,1]:= A1: Display[2,2]:= A2: Display[2,3]:= A3:
Display;

    

If A and B are vectors, matrices, arrays, etc., of the same size and shape, then A *~ B performs the elementwise multiplication.

How would you recommend that this feature be indexed in the help pages, since you seem to have tried to find it? It is listed under ?elementwise

In Maple, e^x is exp(x). There is no caret (^) in it.

See ?Student,Calculus1,VolumeOfRevolution and ?Student,Calculus1,VolumeOfRevolutionTutor

First 199 200 201 202 203 204 205 Last Page 201 of 395