Carl Love

Carl Love

28025 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

Your belief about the cause of the error is spot-on correct. It's trivial to take the positive (or negative) square root of both sides of your ODE, so why don't you try that? 

diff(y(x), x) = y(x)^(-3/2*w-1/2)

It can be done like this:

MulNoK:= (k,n)-> local j: mul(a[k]-a[j], j= {$0..k-1, $k+1..n}):
MulNoK(2,4);

      (a[2] - a[0])*(a[2] - a[1])*(a[2] - a[3])*(a[2] - a[4])

Your first error is that IBC must be passed as a set to pdsolve, so pdsolve(..., {IBC}, ...). Your second error is that the numeric value(s) of phi must be specified before the pdsolve call, so pdsolve(..., eval({IBC}, phi= ...), ...). There are some more errors, but that's what I've worked out so far.

Suggestion: Use more spaces in your code so that it's more readable by a human!

In cylindrical coordinates, the cone is z = r and the paraboloid is z = 12 - r^2. We need the intersection of these surfaces, so

solve(r = 12 - r^2);
                 -4, 3

Note that z = sqrt(x^2 + y^2) implies z >= 0, so we only want the r = 3 solution. So, being "inside" the cone means 0 <= r <= 3, or in Maplese, r= 0..3.

This leads to the plot commands

plots:-display(
    plot3d( #the cone
        [r, theta, r], theta= -Pi..Pi, r= 0..11,
        coords= cylindrical, transparency= 0.7
    ),
    plot3d( #the paraboloid
       [r, theta, 12-r^2], theta= -Pi..Pi, r= 0..3, 
       coords= cylindrical
    )
);

To compute the surface area, I also use cylindrical coordinates: 

f:= (x,y)-> 12-x^2-y^2:
dS:= sqrt(1+D[1](f)(x,y)^2+D[2](f)(x,y)^2);
SA:= Int(sqrt(1+4*r^2)*r, [r= 0..3, theta= -Pi..Pi]);
value(SA);
      Pi/6*(37^(3/2) - 1)

Note the extra factor of r in the integral due to cylindrical coordinates.

How you would "create" the group depends largely on what you want to do with the group. To some extent, all that's needed is the set of elements and the group operation defined on G. Below, I've used &* as the group operation, and I show that the multiplicative group modulo 42 is not cyclic[*1] but that it is a direct product of two cyclic subgroups. In other words, it can be generated by 2 elements, but not by 1. 

restart:
G:= select(x-> igcd(x,42)=1, {$1..41});
      G := {1, 5, 11, 13, 17, 19, 23, 25, 29, 31, 37, 41}
`&*`:= (x,y)-> x*y mod 42:
`&^`:= (x,k::integer)-> x^k mod 42:
ord:= proc(x) local k; for k do until x&^k=1; k end proc:  
ord~([G[]]);
              [1, 6, 6, 2, 6, 6, 6, 3, 2, 6, 3, 2]
Show that G can be generated by two elements:
{seq(seq((5&^i)&*(13&^j), i= 1..6), j= 1..2)};
         {1, 5, 11, 13, 17, 19, 23, 25, 29, 31, 37, 41}
But I can't replace 13 with just any element of order 2:
{seq(seq((5&^i)&*(41&^j), i= 1..6), j= 1..2)};
                     {1, 5, 17, 25, 37, 41}
The reason that it doesn't work is that 41 is a power of 5 (mod 42): 5^3 mod 42 = 41. 

I've made the computational techniques above very simplistic in order to help your understanding. These techniques would be very inefficient for a large modulus.

[*1] It is well known that the multiplicative group modulo n is cyclic iff = 1, 2, 4, pk, or 2*pk for p odd prime.

Obviously the OP Rodey Genius is not interested in learning. However, I had this Answer already 80% done at the time that they deleted the Question, so I'd rather not waste what I coded. Maybe someone else could learn from it.

I chose to code determinants by cofactor expansion because it can be used to illustrate a great many features of Maple programming in one short procedure. The algorithm is straightforward to code recursively. I optimized it the way one would do by hand for a sparse matrix: At each iteration, I expand along the row or column with the most 0s.

I think that even most Maple experts will learn something by studying this procedure.

As a bonus, I found an (admittedly, rather obscure) class of matrices for which my procedure significantly outperforms default methods: half-sparse matrices of bivariate rational functions.

Determinant:= proc(
    A::Matrix(square), 
    {method::identical(cofactors, other):= 'other'}
)
uses AT= ArrayTools, St= Statistics, LA= LinearAlgebra;
local k, n, r, c, T, transp, rm1, rp1, a;
    if method='other' then return LA:-Determinant(A) fi;
    
    # method = cofactors:
    n:= upperbound(A)[1];
    if n<2 then return `if`(n=0, 1, A[1,1]) fi; #base cases
    #Count number of non-zero entries in every row and column:
    T:= St:-Tally~([((r,c):= AT:-SearchArray(A))]);
    if nops~(T) <> [n,n] then return 0 fi; #There's a row or column of 0s.
    T:= sort~(T, 'key'= rhs);
    transp:= evalb(rhs(T[1][1]) > rhs(T[2][1])); #expand a row, or a column?
    r:= lhs(T[`if`(transp, 2, 1)][1]);  rm1:= r-1;  rp1:= r+1;
    add(
        if (a:= `if`(transp, A[k,r], A[r,k])) = 0 then 0
        else
            (-1)^(r+k)*a*thisproc(
                A[[[..rm1,rp1..],[..k-1,k+1..]][`if`(transp, [2,1], [1,2])][]], 
                _options
            )
        fi, 
        k= 1..n
    )
end proc
:
#Test case: half-sparse matrix of bivariate rational functions.
A:= LinearAlgebra:-RandomMatrix(
   9, density= 0.5, 
   generator= (()-> randpoly(x, degree=1)/randpoly([x,y], degree= 1))
):
d1:= CodeTools:-Usage(Determinant(A, method= other)):
memory used=23.74MiB, alloc change=0 bytes, cpu time=938.00ms,
real time=482.00ms, gc time=0ns
d2:= CodeTools:-Usage(Determinant(A, method= cofactors)):
memory used=11.99MiB, alloc change=0 bytes, cpu time=109.00ms, 
real time=130.00ms, gc time=0ns
expand(numer(d1-d2)); #Verify equality of results.
                               0
#Verify results are not identically 0:
eval(d2, [x= 0, y= 0]);
           16003080668688229738162772913864379992367
           -----------------------------------------
            86804823806332323196340023957536768000  

 

There's nothing wrong with Kitonum's Answer, but this alternative form is often useful, especially when n is not integer: For xy, and n real,

(x+I*y)^n = (x^2+y^2)^(n/2)*exp(I*n*arctan(y,x))

Easy analytic geometry solution:

A:= <1,2,3>:  v:= <1,2,2>:  C:= <10,11,12>:
L:= A+t*v:

We want the point where line L intersects the plane normal to containing point C:

B:= eval(L, t= solve(L.v = C.v));

Plot:

plots:-display(
    plots:-implicitplot3d(
        <x,y,z>.v = C.v, ([x,y,z]=~ 0..15)[], 
        style= surface, transparency= .2
    ), 
    plots:-pointplot3d([A,B,C], connect, color= red, thickness= 3),
    orientation= [80, 80, 130], axes= frame, labels= [``$3]
);


 

NULL

n:= 3:

V:= 'LinearAlgebra:-RandomVector(n)' $ n;

V := Vector(3, {(1) = 92, (2) = -31, (3) = 67}), Vector(3, {(1) = 99, (2) = 29, (3) = 44}), Vector(3, {(1) = 27, (2) = 8, (3) = 69})

A:= `<|>`(V); # `<|>` is the columnwise juxtaposition operator.

Matrix(3, 3, {(1, 1) = 92, (1, 2) = 99, (1, 3) = 27, (2, 1) = -31, (2, 2) = 29, (2, 3) = 8, (3, 1) = 67, (3, 2) = 44, (3, 3) = 69})

R:= A^%T.A; # ^%T is the transposition operator.

Matrix(3, 3, {(1, 1) = 13914, (1, 2) = 11157, (1, 3) = 6859, (2, 1) = 11157, (2, 2) = 12578, (2, 3) = 5941, (3, 1) = 6859, (3, 2) = 5941, (3, 3) = 5554})

X:= 1/A; #inverse matrix

Matrix(3, 3, {(1, 1) = 1649/327244, (1, 2) = -5643/327244, (1, 3) = 9/327244, (2, 1) = 2675/327244, (2, 2) = 4539/327244, (2, 3) = -1573/327244, (3, 1) = -3307/327244, (3, 2) = 2585/327244, (3, 3) = 5737/327244})

#Verification that X has the property you asked for:
seq(A.X[..,k], k= 1..3);

Vector(3, {(1) = 1, (2) = 0, (3) = 0}), Vector(3, {(1) = 0, (2) = 1, (3) = 0}), Vector(3, {(1) = 0, (2) = 0, (3) = 1})

``


 

Download BasicLA.mw

If you use option matrixview, then the vertical-axis tickmarks are indexed -1, -2, ..., -20.

The command SPolynomial is in package Groebner not package Ore_algebra.

 

You need to do two things:

1. Your plot range 0..5 is much too big for this problem. Change it to approximately 0..0.01 (for all plots).

2. Add the option maxfun= 0 to the dsolve command. This will remove any limitation on the number of points that it needs to compute.

3. (Optional) Remove output= listprocedure. It doesn't do anything beneficial in this case, and I'm always suspicious of it.

Here it is. Obviously the placement of the sliders could be prettier, but I think mathematically it's all done. MaplePrimes can't handle Explore plots, but it should be there when you download the worksheet.
 

S-I-R infectious disease model with added vaccination-rate parameter

restart
:

Digits:= 15
:

Both I and gamma are difficult to use in Maple. I've changed these to H and y

ODEs:=.
    (diff(S(t),t) = (S+R)(t)*b - S(t)*(beta*H(t) + v + d)),
    (diff(H(t),t) = H(t)*(beta*S(t) - delta - y + b)),
    (diff(R(t),t) = v*S(t) - R(t)*d + y*H(t))
:

ICs:= S(0) = S__0, H(0) = I__0, R(0) = R__0:

SIR:= dsolve(
   {ODEs, ICs}, parameters= [b, beta, v, d, delta, y, S__0, I__0, R__0],
   numeric, maxfun= 0, abserr= 1e-13
):

Explore(
    proc()
    local
        #population growth function scaled to initial pop. = 1 = 100%:
        N:= t-> exp((b-delta-d)*t)
    ;
        SIR('parameters'=
                [b, beta, v, d, delta, y, S__0, I__0, R__0]
        );
        plots:-odeplot(
            SIR, `[]`~(t, [S, H, R](t)*~exp((b-d-delta)*t)), t= 0..50,
            color= [green, yellow, pink], legend= [` S `,` I `,` R `],
            labels= [:-t, `Pop.`], view= [DEFAULT, 0..1]
        )
    end proc(),
    'parameters'= [
        S__0= 0.0..1.0, I__0= 0.0..1.0, R__0= 0.0..1.0,
        b=  0.0..1.0, delta= 0.0..1.0, d= 0.0..1.0, v= 0.0..1.0,
        y= 0.0..0.8, beta= 0.0..2.0
    ]
);
            

 


 

Download VaccinationRateSIR.mw

You're missing a closing parenthesis at the end of the line immediately above end. Note that when your cursor is immediately after a closing bracket, then a faint box appears around it's opening match. 

Looks like you're doing the elliptic curve addition! Good.

You have a few more errors that will just cause weirdness rather than an error message:

  1. Assignment is := not =.
  2. Expand does nothing without mod.
  3. Multiplication should use an *.

Here is a phase portrait for the three-tree-frog system using some parameter values and initial conditions that I just made up. You'll need to do a lot of playing around with the parameters to figure out the rest of part (c) of the problem.
 

restart
:

Digits:= 15
:

From: Steven Strogatz, Nonlinear Dynamics and Chaos, page 304, problem 8.6.9.

 

I do part a) in my head, rewriting the equations in terms of phi and psi:

ODEs:=
   diff(phi(t),t) = -2*H(phi(t))-H(phi(t)+psi(t))+H(psi(t)),
   diff(psi(t),t) = H(phi(t))-H(phi(t)+psi(t))-2*H(psi(t))
;

diff(phi(t), t) = -2*H(phi(t))-H(phi(t)+psi(t))+H(psi(t)), diff(psi(t), t) = H(phi(t))-H(phi(t)+psi(t))-2*H(psi(t))

Part c) The interaction function:

H:= x-> a*sin(x)+b*sin(2*x):

Make up some initial conditions:

ICs:= phi(0) = c, psi(0) = d:

Sol:= dsolve(
   {ODEs, ICs}, parameters= [a,b,c,d],
   numeric, maxfun= 0, abserr=1e-13
):

Experiment with some parameter values:

Sol(parameters= [a=2,b=5,c=2,d=3]):

Phase portrait:

plots:-odeplot(
   Sol, [phi(t),psi(t)], t= 0..50,
   numpoints= 10000, gridlines= false
);

 


 

Download ThreeTreeFrogs.mw

First 121 122 123 124 125 126 127 Last Page 123 of 395