vv

5479 Reputation

9 Badges

3 years, 179 days

MaplePrimes Activity


These are answers submitted by vv

Maple can use theoretically up to 3.8*10^10 digits - if you have enough memory :-).
Your number has 83196715617456  digits (i.e. > 8*10^13) , so it will "overflow".

Using Maple + some elementary math yo can compute

a) the first n digits
b) the last n digits
for any reasonable n.
Try it.

A much more dificult task (probably impossible) would be to compute e.g. the 2 digits in the middle for this number.
If you can do it, please let us know!

 

1.  Use \\  or  /  as  directory separator.

2. The assignment operator is := 
so, use its := 1;

simplify((1/2)!);
  

restart;

# So, u = u(y,t), the rest are constants

pde := diff(u(y, t), t)-1/2*(diff(u(y, t), y, y))-A*sin(M*x-t) = 0;
bc[1] := u(0, t) = 0;
bc[2] := u(10, t) = A*sin(M*x-t);
sys := [pde, bc[1], bc[2]];

diff(u(y, t), t)-(1/2)*(diff(diff(u(y, t), y), y))-A*sin(M*x-t) = 0

 

u(0, t) = 0

 

u(10, t) = A*sin(M*x-t)

 

[diff(u(y, t), t)-(1/2)*(diff(diff(u(y, t), y), y))-A*sin(M*x-t) = 0, u(0, t) = 0, u(10, t) = A*sin(M*x-t)]

(1)

sol:=pdsolve(sys);

u(y, t) = (1/2)*(sin(10+10*I)*(((-1+I)+cos(10-10*I))*sin((1-I)*y)-sin(10-10*I)*cos((1-I)*y))*exp(-I*(M*x-t))-(((1+I-cos(10+10*I))*sin((1+I)*y)+cos((1+I)*y)*sin(10+10*I))*exp(I*(M*x-t))-2*cos(M*x-t)*sin(10+10*I))*sin(10-10*I))*A/(sin(10+10*I)*sin(10-10*I))

(2)

pdetest(sol,sys);

[0, 0, 0]

(3)

You cannot simplify to 0 because it is not 0. But it is a constant.
You can check this in Maple, computing

diff(expr, x);   # 0

where expr is your expression. Do not forget to use u(x) and v(x)  instead of u, v.
 

Tunnell:=proc(N::posint)
local A, a,t, n, isq;
n:=mul(map(`^`~ @ `op`, map(t->[t[1],irem(t[2],2)], ifactors(N)[2])));
isq:=proc(n) local r:=isqrt(trunc(n)); `if`(r^2>n,r-1,r) end;

A:=proc(a,b,c)
  local x,y,z, r:=0;
  for z from 0 to isq(n/c) do
  for y from 0 to isq((n - c*z^2)/b) do
    x:=(n-c*z^2-b*y^2)/a;
    if issqr(x) then r:=r+`if`(z=0,1,2)*`if`(y=0,1,2)*`if`(x=0,1,2) fi
  od od: 
r 
end;

if    n::odd  and 2*A(1,2,32)=A(1,2,8) 
   or n::even and 2*A(2,8,64)=A(2,8,16)
then true else false fi
end:

select(Tunnell, [seq(1 .. 200)]);

[5, 6, 7, 13, 14, 15, 20, 21, 22, 23, 24, 28, 29, 30, 31, 34, 37, 38, 39, 41, 45, 46, 47, 52, 53, 54, 55, 56, 60, 61, 62, 63, 65, 69, 70, 71, 77, 78, 79, 80, 84, 85, 86, 87, 88, 92, 93, 94, 95, 96, 101, 102, 103, 109, 110, 111, 112, 116, 117, 118, 119, 120, 124, 125, 126, 127, 133, 134, 135, 136, 137, 138, 141, 142, 143, 145, 148, 149, 150, 151, 152, 154, 156, 157, 158, 159, 161, 164, 165, 166, 167, 173, 174, 175, 180, 181, 182, 183, 184, 188, 189, 190, 191, 194, 197, 198, 199]

Edit. The proc was optimized for speed.

A module is very convenient and has advantages, but it is not mandatory here.

with(Statistics):

Levy:= (m,s) -> Distribution(
                      PDF=(x->piecewise(x<0, 0, sqrt(s/2/Pi)*exp(-s/2/(x-m))/(x-m)^(3/2))), 
                      CDF= (x->erfc(sqrt(s/2/(x-m)))), 
                      RandomSample = (N->m +~ s /~ ( Quantile~(Normal(0, 1), Sample(Uniform(1/2, 1), N)) )^~2)
                    ):

X := RandomVariable(Levy(m,s)):
PDF(X, x);
Sample(X, 5);
Y := RandomVariable(Levy(1,2)):
PDF(Y, x);
Sample(Y, 6);

EllipticF(z, k)  has a bug for  z --> infinity  for any k (not necessarily complex).

It seems to be an old and known one.

Let's illustrate for k=1/2.

 

f:= EllipticF(z,1/2);

EllipticF(z, 1/2)

(1)

limit(f, z=infinity); # fails

limit(EllipticF(z, 1/2), z = infinity)

(2)

MultiSeries:-limit(f, z=infinity); # wrong

4*EllipticK(2)-(2*I)*EllipticK((1/2)*(-3)^(1/2)*4^(1/2))

(3)

evalf(%);

3.371500710-6.469546944*I

(4)

eval(f, z=1e6);

0.2000000000e-5-2.156515648*I

(5)

# The correct limit is

- EllipticCK(1/2)*I;
evalf(%);

-I*EllipticCK(1/2)

 

-2.156515647*I

(6)

# asympt is also wrong

asympt(f, z, 2); # fails

Error, (in asympt) unable to compute series

 

MultiSeries:-asympt(f, z, 2); # wrong

4*EllipticK(2)-(2*I)*EllipticK((1/2)*(-3)^(1/2)*4^(1/2))-2/z+O(1/z^3)

(7)

evalf(%);

3.371500710-6.469546944*I-2./z+O(1/z^3)

(8)

restart;

with(GraphTheory):

A0:=<
0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 1, 1, 0, 1, 0, 0, 0, 1;
1, 1, 0, 0, 1, 0, 1, 1, 1, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0, 0
>;

Matrix(4, 10, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 1, (2, 5) = 0, (2, 6) = 1, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 1, (3, 1) = 1, (3, 2) = 1, (3, 3) = 0, (3, 4) = 0, (3, 5) = 1, (3, 6) = 0, (3, 7) = 1, (3, 8) = 1, (3, 9) = 1, (3, 10) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0})

(1)

A:=Matrix(10,A0);

Matrix(10, 10, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 1, (2, 5) = 0, (2, 6) = 1, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 1, (3, 1) = 1, (3, 2) = 1, (3, 3) = 0, (3, 4) = 0, (3, 5) = 1, (3, 6) = 0, (3, 7) = 1, (3, 8) = 1, (3, 9) = 1, (3, 10) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (5, 7) = 0, (5, 8) = 0, (5, 9) = 0, (5, 10) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0, (6, 7) = 0, (6, 8) = 0, (6, 9) = 0, (6, 10) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = 0, (7, 6) = 0, (7, 7) = 0, (7, 8) = 0, (7, 9) = 0, (7, 10) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = 0, (8, 6) = 0, (8, 7) = 0, (8, 8) = 0, (8, 9) = 0, (8, 10) = 0, (9, 1) = 0, (9, 2) = 0, (9, 3) = 0, (9, 4) = 0, (9, 5) = 0, (9, 6) = 0, (9, 7) = 0, (9, 8) = 0, (9, 9) = 0, (9, 10) = 0, (10, 1) = 0, (10, 2) = 0, (10, 3) = 0, (10, 4) = 0, (10, 5) = 0, (10, 6) = 0, (10, 7) = 0, (10, 8) = 0, (10, 9) = 0, (10, 10) = 0})

(2)

G:=Graph(A);  # directed graph

GRAPHLN(directed, unweighted, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], Array(%id = 18446744074365436142), `GRAPHLN/table/1`, 0)

(3)

DrawGraph(G, style=circle);

 

node:=2;
A||node:=Matrix(10):  A||node[node]:=A[node]:
G||node:=Graph(A||node):
DrawGraph(G||node, style=circle);

2

 

 


Download DirGraph.mw

Why do you think that the desired solution must be obtained for _C1 = 0?
Use e.g.
S:=dsolve({E1, phi(0) = b*a/(a*c-b^2)},  phi(xi));
symplify([S], symbolic);

You will find your solution, but be careful about branches.

(Note also that the ODE is not in an explicite form, so even with an initial condition it could have infinitely many solutions).

restart;

xnextHelley:= x - 2*f(x)*D(f)(x)/( 2* D(f)(x)^2 - f(x)*D(D(f))(x) ); # the method

x-2*f(x)*(D(f))(x)/(2*(D(f))(x)^2-f(x)*((D@@2)(f))(x))

(1)

f:= x -> a*x + b*x^2 + c*x^3 + d(x)*x^4;  # WLOG, x[infinity]=0

proc (x) options operator, arrow; a*x+b*x^2+c*x^3+d(x)*x^4 end proc

(2)

series(xnextHelley, x);

series(-((a*c-b^2)/a^2)*x^3-((3*a*d(0)-3*b*c-3*(a*c-b^2)*b/a)/a^2)*x^4-((6*a*(D(d))(0)-3*b*d(0)-3*c^2-(1/2)*(a*c-b^2)*(6*a*c+6*b^2)/a^2-9*(a^2*d(0)-2*b*c*a+b^3)*b/a^2)/a^2)*x^5+O(x^6),x,6)

(3)

# So,

abs(x[n+1] - x[infinity]) <= M * abs(x[n] - x[infinity])^3;  # ==>  order=3

abs(-x[n+1]+x[infinity]) <= M*abs(x[n]-x[infinity])^3

(4)

 

Note first that A^n can be computed symbolically using LinearAlgebra:-MatrixPower(A,n);
The entries are the Fibonacci numbers, and as Rouben said, it is easy to see that the maximum entry is fibonacci(n+1) at (1,1).
Here is a procedure to compute the inverse Fibonacci numbers.

inversefib:=proc(y::posint)
local n1,n2, s1,s2;
if y=1 then return 'Exact'=2 fi;
n1:=ln( (1/2)*sqrt(5)*y+(1/2)*sqrt(5*y^2+4) )/ln((sqrt(5)+1)/2); # n even
n2:=ln( (1/2)*sqrt(5)*y+(1/2)*sqrt(5*y^2-4) )/ln((sqrt(5)+1)/2); # n odd
s1,s2:=simplify~([n1,n2])[];
if   s1::even then  'Exact'=s1
elif s2::odd  then  'Exact'=s2
else  'Next'= ceil(s1) fi
end proc:

Examples.
inversefib(2019);
    
Next = 18
inversefib(3524578);
     Exact = 33
inversefib(10^12);
     Next = 60

 

 

remove[flatten](member, A, B);

If you want it inplace, use

remove[flatten,inplace](member, A, B);

plot(arctan(x),x=-20..20, 
title=typeset(
H__0,": ", sigma[1]^2=sigma[2]^2, "     ",  
H__1,": ", sigma[1]^2<>sigma[2]^2, "\n", MyTitle)
);

Yes, it seems that

seq(C[i+1]-C[i], i = 1 .. N-1):
does not finish the job.
But,
seq(C[i+1]-i^2, i = 1 .. N-1):
works.
Probably a bug.

 

 

1 2 3 4 5 6 7 Last Page 1 of 64