Carl Love

Carl Love

28015 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

See the command Finance:-BrownianMotion (help page ?Finance,BrownianMotion).

There are three major problems with your code:

1. In the piecewise command, the boolean conditions (i.e., the inequalities in this case) must come before their corresponding expressions, even though they appear after in the standard forms of display (whether in Maple or just typeset mathematics).

2. You can't use x both as a vector and as a symbolic variable (such as a variable of integration). Below, I've changed the symbolic x to xi.

3. Variables that appear in expressions after an arrow -> are not evaluated at the time the function (procedure, arrow expression) is defined. They aren't evaluated until the function is called. This could be corrected with unapply (as suggested by VV above), but I think that it's better to eliminate the loop entirely.

To correct these problems, delete the entire loop where phi is defined. Replace it with this single function definition:

phi:= (k::posint, X::Vector(realcons), h::algebraic)->
    x-> piecewise(Or(x < X[k], x > X[k+2]), 0, x <= X[k+1], x-X[k], X[k+2]-x)/h
:

And change the int command to

F[k]:= int(sin(Pi*xi)*phi(k,x,h)(xi), xi= 0..1)

It's quite simple, and nothing needs to be done implicitly:

ode:= diff(f(x),x)*exp(f(x)-x^2-1) = 2*x:
dsolve({ode, f(0)=1});
                                 2    
                         f(x) = x  + 1
int(rhs(%), x= 0..1);
                               4
                               -
                               3

By the way, note that the ODE is separable as 

int(exp(y), y) = C+int(2*x*exp(x^2+1), x)

where y= f(x).

The Optimization package is for numeric solutions. It can't give solutions for one variable in terms of other variables.

The radius for 4 points is obviously 1/sqrt(2). A set of four points satisfying the conditions is {[1/sqrt(2), 0], [0, 1/sqrt(2)], [-1/sqrt(2), 0], [0, -1/sqrt(2)]}.

For any n, the radius r can be determined by solving the equation

sol:= solve(1=(r-r*cos(2*Pi/n))^2+(r*sin(2*Pi/n))^2, r)

and taking the positive branch of the square root. The result (very easy to obtain by hand also) is

r = 1/sqrt(2 - 2*cos(2*Pi/n))

That can be simplified via a half-angle identity to 

r = csc(Pi/n)/2

A list of n concyclic points satisfying the condition is then

[seq](r*~[cos,sin](2*Pi*k/n), k= 1..n) 

with as given above.

Any three non-collinear points uniquely determine a circle (and, of course, a triangle). (That existence and uniqueness of the circle was surprising to me when I first learned of it, so many decades ago.) So, we just need the center and radius of the circle determined by the vertices of the triangle. This can be done easily with the geometry package. Below, I generate three points at random, find the circle, and plot the circle, triangle, and circumcenter.

restart:
R:= rand(-9..9):
G:= geometry:
G:-circle(C1, G:-point~([A,B,C], '['R()'$3]'$2)):
Triang:= G:-coordinates~([A, B, C, A]);
         Triang := [[5, -8], [3, 2], [-4, -6], [5, -8]]

Ccen:= G:-coordinates(G:-center(C1));
                              [109  -305]
                      Ccen := [---, ----]
                              [86    86 ]
R:= G:-radius(C1);
                      1         (1/2)     (1/2)
                R := ---- 124865      3698     
                     3698                      
plot(
    [[(R*~[cos,sin](t) +~ Ccen)[], t= -Pi..Pi], Triang, [Ccen]], 
    style= [line$2, point], symbol= solidcircle, scaling= constrained,
    thickness= 5
);

If you want to work through (with Maple) the algebra that underlies the computations done by geometry above, let me know; it's not terribly complicated.

One way is like this:

x:= <1, 2, 3, 4>;
k:= Matrix(4, 4, (i,j)-> x[i]+x[j]);

Nested for loops could also be used, but it's not ideal.

Do you want the matrix to contain only the coefficients of the power-3 terms of the equations? In that case, do

KKNL:= Matrix((nops(EqML), nops(Var[4])), (i,j)-> coeff(EqML[i], Var[4][j]))

Acer's Answer is good. I just want to add what I think may have contributed to your original confusion about type indexed versus type symbol: When option symbol= a is used with a Matrix command, the nonzero entries of the resulting matrix are type indexed, not type symbol. It is the a (alone) which is the symbol; not the entries such as a[1,1].

Your desired form of the solution contains several factors that don't depend on t. These factors could be absorbed into the constants of integration without changing the correctness of the solution.

I agree with @tomleslie that you should use add instead of sum and local indices and shouldn't use unevaluation quotes. But I think the more-fundamental problem is using catenation rather than true indexing (such as A[i,j,k]). So, I'd do this (in 1D input):

p:= (i,j)-> local k; add(A[i,j,k], k= 1..3):
q:= i-> local j; add(p(i,j), j= 1..3):
r:= add(q(i), i= 1..3);

I can import and parse your 2800x2800 matrix of rational functions that you posted as file "f.txt" in only 7 seconds, like this:

restart:
M:= CodeTools:-Usage(
    subsindets(
        ImportMatrix("C:/Users/carlj/Downloads/f.txt"),
        string,
        parse
    )
):
memory used=185.79MiB, alloc change=195.61MiB, 
cpu time=6.69s, real time=6.18s, gc time=1.22s

upperbound(M);
                           2800, 2800

#Example entry:
M[9,9];
          /                                         
    1     |      2                                  
--------- |a[1] h  (K__ref Kuy0[1] + K__ref Kuyb[1])
 3        \                                         
h  K__ref                                           

                      3                                            
     144 a[1] K__ref h  #mover(mi("A"),mo("&uminus0;"))[6, 6, 1]   
   + ----------------------------------------------------------- - 
                                b[1]                               

        1       /     
  ------------- \a[1] 
  17 b[1] ro[1]       

                                                            2 
  #mover(mi("&omega;",fontstyle = "normal"),mo("&uminus0;"))  

          2        \
  K__ref h  I__0[1]/

                                                 \
     1        2                                  |
   + -- b[1] h  (K__ref Kux0[1] + K__ref Kuxa[1])|
     17                                          /


The MathML operators shown above (movermomi, etc.) won't appear in ordinary prettyprinted output. Their appearence above is an artifact of posting on MaplePrimes.

The subsindets(..., string, parse) is needed because symbolic expressions are being imported as strings, which isn't really a problem because parse very quickly converts them back to symbolic expressions.

I "cleaned up" your code, as you suggested.

You mentioned the possibility of using ithprime. Even better is nextprime. It's not really an optimization, but it leads to cleaner code.

Your idea of prime pairs with a fixed difference can be easily extended to prime constellations (aka tuples). For example, [11, 13, 17, 19] is a 4-tuple with difference pattern [2, 6, 8].

restart:
PrimeTuples:= proc(
    diffs::list(posint),
    #Three keyword-option stop criteria:
    {
        Interval::range({posint, infinity}):= 2..infinity,
        Timelimit::positive:= infinity,
        Count::{posint, infinity}:= infinity
    }
)
local
    P:= rtable(1..0), k:= 0, q:= lhs(Interval)-1,

    #Append next tuple to return array:
    Tuple:= proc()
    local T;
        q:= nextprime(q); T:= q +~ diffs;
        if andmap(isprime, T) then k:= k+1; P(k):= [q, T[]] fi;
        return
    end proc,

    #Check for cases with provably finitely many tuples: 
    Finite:= proc()
    local p:= 1, nd:= nops(diffs) + 1, df:= {diffs[]};
        do
            p:= nextprime(p);
            if p > nd then break fi;
            if nops(irem~(df, p) minus {0}) = p-1 then
                while q < p do Tuple() od; return true
            fi
        od;
        false
    end proc,

    #Check for cases with likely infinitely many tuples:        
    Infinite:= proc()
    option procname;
    local upper:= rhs(Interval);
        if {upper, Count, Timelimit} = {infinity} then
            error "no stop criterion"
        fi;
        while q < upper and k < Count do Tuple() od;
        return
    end proc       
;    
    if not Finite() then    
        try timelimit(eval(Timelimit, infinity= -1), Infinite())
        catch "time expired":
        end try
    fi;

    convert(P, list)
end proc
: 

PrimeTuples([2,4]);
                          [[3, 5, 7]]

PrimeTuples([2,8,14,26]);
            [[3, 5, 11, 17, 29], [5, 7, 13, 19, 31]]

PrimeTuples([2,4,6]);
                               []

PrimeTuples([2,6,8]);
Error, (in PrimeTuples) no stop criterion

PrimeTuples([2,6,8], Count= 9);
   [[5, 7, 11, 13], [11, 13, 17, 19], [101, 103, 107, 109], 
    [191, 193, 197, 199], [821, 823, 827, 829], [1481, 1483, 1487, 1489], 
    [1871, 1873, 1877, 1879], [2081, 2083, 2087, 2089], [3251, 3253, 3257, 3259]]

PrimeTuples([48], Interval= 2..1000);
[[5, 53], [11, 59], [13, 61], [19, 67], [23, 71], [31, 79], 
  [41, 89], [53, 101], [59, 107], [61, 109], [79, 127], 
  [83, 131], [89, 137], [101, 149], [103, 151], [109, 157], 
  [131, 179], [149, 197], [151, 199], [163, 211], [179, 227], 
  [181, 229], [191, 239], [193, 241], [223, 271], [229, 277], 
  [233, 281], [263, 311], [269, 317], [283, 331], [311, 359], 
  [331, 379], [349, 397], [353, 401], [373, 421], [383, 431], 
  [401, 449], [409, 457], [419, 467], [431, 479], [439, 487], 
  [443, 491], [461, 509], [499, 547], [509, 557], [521, 569], 
  [523, 571], [569, 617], [571, 619], [593, 641], [599, 647], 
  [613, 661], [643, 691], [653, 701], [661, 709], [691, 739], 
  [709, 757], [739, 787], [761, 809], [773, 821], [809, 857], 
  [811, 859], [829, 877], [839, 887], [859, 907], [863, 911], 
  [881, 929], [919, 967], [929, 977], [971, 1019], [983, 1031], 
  [991, 1039]] 


Download PrimeTuples.mw

There is a stock command for it: Fractals:-EscapeTime:-Newton. See its help page ?Fractals,EscapeTime,Newton.

The 2nd partial derivatives with respect to the 1st independent variable should be D[1,1](W)(0, tau) and D[1,1](W)(inf, tau).

First 53 54 55 56 57 58 59 Last Page 55 of 394