Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

I saw that you just asked a Question about the derivative of Eq. 6. That Question got deleted, I guess because it's too similar to this one.

Anyway, in that Question, you took the derivative of both sides of Eq6 to get Eq6_2, and you asked for a plot of its solution. This operation is redundant or superfluous because taking the derivative of both sides of an ODE doesn't fundamentally change its solution (provided that an appropriate initial condition is also added). I suspect that what you really want is a plot of the derivative of f(eta), which is a different thing than the solution of the derivative of an equation that contains f(eta). To plot the derivative of f(eta), all that you need to do is change f(eta) to diff(f(eta), eta) in the odeplot command. No other change is needed. Doing that, and also removing the plot title, I get

Here is a complete example of a pseudo-random number generator (PRNG), just intended as an educational example so that you can see how they work. It's surprisingly short and simple. It generates integers between 0 and 100, inclusive. The only difference between this and an actual PNRG is that the actual one would use much larger integers and probably something a little bit more complicated than a linear modular congruence. The way that the Randomize in this module works is almost identical to Maple's own randomize. Note that floats are generated from integers. That's the way it actually works in Maple.

Rand:= module()
local 
    A:= 47, B:= 53, M:= 101, #Note that these are prime.
    seed,
    ModuleApply:= proc() seed:= irem(A*seed+B, M) end proc,
    ModuleLoad:= proc() seed:= 1 end proc
;
export 
    Randomize:= proc(s::posint:= 0) seed:= `if`(s=0, Value(Time()), s) end proc,
    Float:= ()-> evalf(Rand()/M)
;
    seed:= 1
end module
:    
'Rand()' $ 9;
                100, 6, 32, 42, 7, 79, 29, 2, 46
Seed1:= Rand:-Randomize();
                     Seed1 := 1577955051735
'Rand()' $ 9;
               34, 35, 82, 69, 64, 31, 96, 20, 84
'Rand()' $ 9;
               62, 38, 21, 30, 49, 33, 89, 95, 74
Rand:-Randomize(Seed1): #Restart sequence from Seed1.
'Rand()' $ 9;
               34, 35, 82, 69, 64, 31, 96, 20, 84

 

These 3 boundary conditions inside the loop are nonsense:

D(F[k])(last) = 0, H[k](last) = 0, T[k](last) = 0

They're nonsense because last is the number of loop iterations; it's not an endpoint of a real interval.

Yes, compiling on each node helps. Here is an example. The part about wasting time is just so that the task is nontrivial and thus I can measure the efficiency of the multiprocessing and know for sure that multiple processors were used.
 

restart:

Sq:= proc(x::integer[4])::integer[4];
local i::integer[4]:= 0, k::integer[4];
    #This loop is only to waste time (approx 4 secs):
    for k to 2\000\000\000 do i:= i+1 od;
    x^2
end proc:

Sq_comp:= Compiler:-Compile(Sq):

#Measure time for one simple run:
rt:= CodeTools:-Usage(Sq_comp(1), output= realtime);

memory used=1.30KiB, alloc change=0 bytes, cpu time=4.34s, real time=4.34s, gc time=0ns

4.342

nc:= kernelopts(numcpus);

8

m:= 3: #can be any small posint

S:= [$1..m*nc]: #input data for distributed processing

#
#Here I show that this naive technique generates an error:
st:= time[real]():
R:= Grid:-Map(Sq_comp, S):
rt2:= time[real]() - st:

Error, (in Grid:-Map) invalid input: op expects 1 or 2 arguments, but received 0

#######################
# Try a different way #
#######################

#Save definitions so that I can restart:
currentdir("C:/Users/carlj/desktop"):
save Sq, rt, nc, m, S, "GridTest.mpl"; #not Sq_comp

#
restart:

#
currentdir("C:/Users/carlj/desktop"):
read "GridTest.mpl":

#
#Correct that error by compiling on each node:
Grid:-Set(Sq);
Grid:-Set(Sq_comp= 'Compiler:-Compile'(Sq));

#
#Attempt distributed processing again:
st:= time[real]():
R:= Grid:-Map(Sq_comp, S):
rt2:= time[real]() - st:

#
#Check accuracy:
evalb(R = S^~2);

true

#
#Check efficiency:
E:= m*rt/rt2;

.6957961647

 


 

Download GridComp.mw

There may also be a way to do it without compiling on each node. I don't know.

Like this:

plots:-textplot(
   [seq(
       [X[k], Y[k], sprintf("%3d", k), color= COLOR(HSV, .85*(k-1)/(N-1), 1, .8)],
       k= 1..N
   )],
   font= ["times", "bold", 9], axes= none, size= [700$2]
);

Functions that don't appear in the ODE system can also be plotted with odeplot. So, you can do

plots:-odeplot(Sol, [P,V](x), x= 0..10)

even though only P(x) appears in the differential equation(s). So, there's no need to change your original ODE(s).

The above command plots one curve, V vs P, not two separate curves.

Just to summarize the essential point, which may be difficult to see amidst the dense numeric code: The parameter of animate is inherently real-valued; it must be somehow converted to an integer before being used as an index to your table SOLNSuy. The alternative to animate is to create an indexable container of static plots and pass it to display(..., insequence).

What you found, unapply, is the most widely recommended way to do it.

The solution is simply the obvious: f(x/a).

I don't understand exactly why what you had doesn't work, but I strongly prefer using eval rather than assign. By simply making that change, it now works. Replace your last line of code with

Sol:= {};
for i from 0 to N do   
    Sol:= Sol union {dsolve({cond[1][i], eval(equ[2][i], Sol)}, {f[i](x)})}
end do;

It can be done like this:

C:= lcoeff(y^2 - x^2/y, order= plex(x,y), 't'):
C*t;

                       -x^2/y

In many top-level Maple commands, the concept of polynomial is generalized a bit to allow for negative exponents. 

You're making two mistakes: The first is mixing up integer and polynomial computation; the second is mixing up exact and floating-point computation. 

Integer vs polynomial: You're making exactly the same mistake as you were making regarding gcd. The mod commands (there's a whole family of them) are intended primarily for polynomials and matrices over finite rings (including finite fields). The command for computing the remainder of integer division is irem. Do you see a pattern here? Top-level commands for pure-integer arithmetic begin with i. Your command x mod 3 returns x because x is a polynomial; it's not a placeholder for an integer.

Exact vs floating-point: Neither of these exact commands is appropriate for use as the first argument to plot, which is necessarily based on floating-point computation. 

Replace sum with add. The problem with sum is that the derivative of u[i] is evaluated before has a value. Since u[i] hasn't been defined in terms of x for generic i, its derivative is 0. On the other hand, add evaluates its index variable to a value in the specified range before evaluating its expression. 

One is tricked into believing that the first assembly in VV's Answer is a right triangle, making the area 32-1/2. But actually that assembly is a quadrilateral. The slope of the hypotenuse of the green triangle is 2/5, while that of the red triangle is 3/8. The difference of the angles is only about 1.2 degrees, or a relative difference of about 6%. Thus the supposed hypotenuse of the overall triangle is actually 2 sides of a quadrilateral, and its area is actually 32.

This is not a bug. The 2 represents the precision at which the computations are done, not the precision of the final result. The number 1014.1 can't be represented adequately in 2 digits.

First 119 120 121 122 123 124 125 Last Page 121 of 395