acer

LUDecomposition...

Answered: acer 32405

November 16 2023

3 1

Is this a correct import of your Matrix?

>

restart;

>	kernelopts(version);

>	with(LinearAlgebra):

>	M := subsindets(ImportMatrix(cat(kernelopts(homedir), "/mapleprimes/VarchenkoMatrix.csv")), string,parse);

_rtable[36893628039329622484]

>	ans := CodeTools:-Usage( simplify(Determinant(LUDecomposition(M,output=U))) );

memory used=9.59MiB, alloc change=0 bytes, cpu time=73.00ms, real time=74.00ms, gc time=0ns

>	f := rand(-37..37): G := [seq(x\|\|i=f(),i=1..5)]: Determinant(eval(M,G)) - eval(ans,G);

>	I find these interesting. andmap(type,LUDecomposition(M,output=U),polynom); add(length~(LUDecomposition(M,output=U))); Normalizer := simplify@factor: Determinant(LUDecomposition(M,output=U)); simplify(%);

-(-x2^2+1)*(-x1^2+1)*(-x4^2+1)*(x3^2-1)^15*(x2^2-1)^15*(x1^2-1)^15*(x4^2-1)^15*(x5^2-1)^15*(-x5^2+1)*(-x3^2+1)

Download det_fun.mw

common "size", or common range&scale&siz...

Answered: acer 32405

November 14 2023

1 5

You can force each to render with a common size (dimensions of the plot window), or you could also have each be displayed with a common range.

Apologies that this site renders such array-plots with an ugly thick black border, etc. That looks better in actual Maple. It's just to show them side-by-side here, and not part of the answer.

>

restart;

>	randomize():

>	with(plottools,getdata): with(plots,display):

>	plots:-setoptions(size=[300,300], gridlines=false):

>

P0 := proc(n, epsilon)
local m, N, D1, Modu, fs, f1, f2, f, q, s, r, t, eps, pp, k;
uses LinearAlgebra, plots;
eps := epsilon;
m := n;
N := eval(RandomMatrix(m, m));
N := N/evalf(Norm(evalf(N), 2));
D1 := DiagonalMatrix(Eigenvalues(N));
pp := pointplot({seq([Re(D1[k, k]), Im(D1[k, k])], k = 1 .. m)});
f := (i, j) -> if i = j then D1[i, i]; elif i = j + 1 then (1 - abs(D1[i, i])^2)^(1/2)*(1 - abs(D1[j, j])^2)^(1/2); elif j + 1 < i then (1 - abs(D1[i, i])^2)^(1/2)*(1 - abs(D1[j, j])^2)^(1/2)*mul(-conjugate(D1[t, t]), t = j + 1 .. i - 1); else 0; end if;
Modu := Matrix(m, (i, j) -> f(i, j)); fs := (x, y) -> Norm(1/((x + y*I)*IdentityMatrix(m) - N), 2) - 1/eps;
f1 := (x, y) -> Norm(1/((x + y*I)*IdentityMatrix(m) - D1), 2) - 1/eps;
f2 := (x, y) -> Norm(1/((x + y*I)*IdentityMatrix(m) - Modu), 2) - 1/eps;
s := implicitplot(fs, -2 .. 2, -2 .. 2, axes = boxed, color = blue, gridrefine = 3, scaling = constrained, resolution = 1000);
q := implicitplot(f1, -2 .. 2, -2 .. 2, axes = boxed, gridrefine = 3, scaling = constrained, resolution = 1000);
r := implicitplot(f2, -2 .. 2, -2 .. 2, axes = boxed, color = green, gridrefine = 3, scaling = constrained, resolution = 1000);
return display({q, s, r, pp});
end proc:

>	A1 := P0(6, 1e-5): A2 := P0(6, 1e-5): A3 := P0(6, 1e-5):

>	A4 := P0(6, 1e-5): A5 := P0(6, 1e-5): A6 := P0(6, 1e-5):

>

# These all have the same "SIZES" (OP's wording).
# Personally I find it hard to interpret and compare
# the results, since the ranges are all different.

opts := size=[300,300], scaling=unconstrained:

display(Array([
[display(A1, opts), display(A2, opts), display(A3, opts)],
[display(A4, opts), display(A5, opts), display(A6, opts)]
]));

>

# Default of this proc is to enlarge by 10%,
# ie. scale the ranges by a factor of 1.1

getR := proc(P::seq(specfunc(anything,PLOT)),
             {scale::And(float,positive):=1.1})
  local xr,yr,xlo,xhi,ylo,yhi,xmed,ymed,xd,yd;
  xr,yr := getdata(plots:-display(P),'rangesonly')[];
  xlo,xhi,ylo,yhi := op(xr),op(yr);
  xmed,ymed := (xhi+xlo)/2,(yhi+ylo)/2;
  xd,yd := xhi-xlo,yhi-ylo;
  xlo,xhi := xmed-scale*xd/2,xmed+scale*xd/2;
  ylo,yhi := ymed-scale*yd/2,ymed+scale*yd/2;
  xlo..xhi, ylo..yhi
end proc:

>	hrng,vrng := getR(A1,A2,A3,A4,A5,A6);

>	opts := view=[hrng,vrng], xtickmarks=evalf((trunc~([0,op(hrng)]~10)/~10)), ytickmarks=evalf((trunc~([0,op(vrng)]~10)/~10));

>	# These all share a common range, which (IMO) # makes them easier to interpret and compare side-by-side display(Array([ [display(A1, opts), display(A2, opts), display(A3, opts)], [display(A4, opts), display(A5, opts), display(A6, opts)] ]));

>

Download rzarouf_01.mw

evala, etc...

Answered: acer 32405

November 13 2023

2 2

Using Maple 2023.1,

expr:=  (x^(-(44 + 12*sqrt(69))^(1/3)/6 + 10/(3*(44 + 12*sqrt(69))^(1/3)) + 2/3)) - (x^(sqrt(69)*2^(1/3)*((11 + 3*sqrt(69))^2)^(1/3)/100 - (11*(44 + 12*sqrt(69))^(2/3))/600 - (44 + 12*sqrt(69))^(1/3)/6 + 2/3)):

evala(expr);

            0

radnormal(expr);

            0

combine(factor(expr));

            0

factor(combine(expr));

            0

fsolve...

Answered: acer 32405

November 13 2023

2 0

f := x-> 2*sin(x) + 1/2*x - 1:

fsolve(f, -2*Pi..2*Pi, maxsols=10);

     0.4090496716, 3.535612202, 5.308993144

example...

Answered: acer 32405

November 12 2023

1 1

>

restart;
n := 0:
for a to 10 do
for b to 10 do
for c to 10 do
for d to 10 do
for t to 10 do
mydelta := 4*a*t^2 - 4*b*d*t + 4*c*d^2 - 4*a*c + b^2;
if 0 < mydelta and type(mydelta, integer) and -d^2 + a <> 0 then
x1 := (2*d*t - b + sqrt(4*a*t^2 - 4*b*d*t + 4*c*d^2 - 4*a*c + b^2))/(2*(-d^2 + a));
x2 := -(-2*d*t + sqrt(4*a*t^2 - 4*b*d*t + 4*c*d^2 - 4*a*c + b^2) + b)/(2*(-d^2 + a));
if 0 < d*x1 + t and 0 < d*x2 + t and type(x1, integer) and type(x2, integer) and nops({a, b, c, d, t}) = 5 and c <> t^2 and -d^2 + a <> 0 and nops({x1, x2}) = 2 then n := n + 1; L[n] := [a, b, c, d, t]; end if; end if; end do; end do; end do; end do;
end do;
L := convert(L, list);
nops(L);

Download minhthien.mw

example...

Answered: acer 32405

November 12 2023

0 1

Your code attempts to evaluate the derivatives of the individual functions as if the result of each solving call were of the same sequence-of-equations form that dsolve (numeric) produces. Eg.

eval( diff(Theta(xi, eta), eta), Ans[j] )

But that's not right -- you are using pdsolve (numeric), and its result (assigned to any Ans[j]) is a module, not a set or list of equations. This is your basic mistake.

You could instead try using the D operator to compute the various derivatives numerically, similarly to what I did here and here, say. (I suspect that you have seen that code, or something very similar to its approach.)

Below I apply the D operator to procedures TTheta, PPhi, and ff which can produce similar results as Ans[k]:-value for Theta(..), Phi(..), and f(..). Below I compare this approach for the Theta plot, as a check. You should check the others.

Note that the numeric computation of the requested derivatives is then computed as numeric approximations, via fdiff as called by evalf(D(...)). I also adjusted a step-size, and you might need to adjust such further. I do not know what accuracy is possible using this approach: using a difference scheme to approximate a derivative of a solution based itself on a finite-difference computation. (At some precision that would likely all break down...)

I did not instead try to augment your DE or figure out more boundary conditions, etc.

Check it over. Do any of the plots appear as you expect?

>

>	for j to numelems(lambdaVals) do Ans[j] := pdsolve((eval([OdeSys, Cond], lambda = lambdaVals[j]))[], numeric, spacestep=0.05, compile=true); end do:

>	TTheta := proc(Xi, Eta, j) if not [Xi, Eta, j]::(list(numeric)) then return ('procname')(args) end if;eval(Theta(xi, eta), (Ans[j]:-value(xi = Xi))(Eta)) end proc:

>	PPhi := proc(Xi, Eta, j) if not [Xi, Eta, j]::(list(numeric)) then return ('procname')(args) end if;eval(Phi(xi, eta), (Ans[j]:-value(xi = Xi))(Eta)) end proc:

>	ff := proc(Xi, Eta, j) if not [Xi, Eta, j]::(list(numeric)) then return ('procname')(args) end if;eval(f(xi, eta), (Ans[j]:-value(xi = Xi))(Eta)) end proc:

>	Ng := proc(xi,eta,j) alphat(1+4Rd(1/3))D[2](TTheta)(xi, eta, j)^2+LtalphacD[2](PPhi)(xi, eta, j)^2/alphat+LtD[2](PPhi)(xi, eta, j)D[2](TTheta)(xi, eta, j)+Br(xiD[2](ff)(xi, eta, j)-sin(xi))^2+(Brxixi)D[2,2](ff)(xi, eta, j)^2-(1/3)BrdeltaBxi^4*D[2,2](ff)(xi, eta, j)^4; end proc:

>

Bj := proc(xi,eta,j) alphat*(1+4*Rd*(1/3))*D[2](TTheta)(xi, eta, j)^2+Lt*alphac*D[2](PPhi)(xi, eta, j)^2/alphat+Lt*D[2](PPhi)(xi, eta, j)*D[2](TTheta)(xi, eta, j)/alphat*(1+4*Rd*(1/3))*D[2](TTheta)(xi, eta, j)^2+Lt*alphac*D[2](PPhi)(xi, eta, j)^2/alphat+Lt*D[2](PPhi)(xi, eta, j)*D[2](TTheta)(xi, eta, j)+Br*(xi*D[2](ff)(xi, eta, j)-sin(xi))^2+(Br*xi*xi)*D[2,2](ff)(xi, eta, j)^2-(1/3)*Br*deltaB*xi^4*D[2,2](ff)(xi, eta, j)^4;
end proc:

>

Download error_in_plot_ac.mw

example...

Answered: acer 32405

November 11 2023

2 4

Is this the kind of thing that you're trying to accomplish?

If it's not quite right then please just mention the detail; there's a decent chance that it could be fixed up.

Download TS_multi_eq_uneval.mw

Note the absence of gray multiplication symbols, ie. they're rendered in the more usual black here.

And you only need to pass in A the assigned name, not all those other variant arguments.

And in the original assignment to A you only need a single pair of wrapping uneval quotes. (Using two pairs made it look like a kludge.)

For fun, with asterisks instead of center-dots for multiplication in the third output term. It's just an example, in case you'd care to play/adjust. It also could likely be done easily with all spaces (denoting multiplication implicitly), or all center-dots throughout all the terms, or all asterisks throughout.

Download TS_multi_eq_uneval_2.mw

add...

Answered: acer 32405

November 09 2023

0 1

I recommend that you use the add command, rather than the sum command, for adding up this finite number of terms.

>	m := < < -0.143 \| -1.145 >, < 0.138 \| -0.832 >, < -0.177 \| -0.149 >, < -0.007 \| -0.028 > >;

Matrix(4, 2, {(1, 1) = -.143, (1, 2) = -1.145, (2, 1) = .138, (2, 2) = -.832, (3, 1) = -.177, (3, 2) = -.149, (4, 1) = -0.7e-2, (4, 2) = -0.28e-1})

>	X:=2: Y:=3:

>	Your mistake, since the Matrix cannot be indexed by symbolic names like i,j (ie. not numbers). This is what would normally get evaluated here, before being passed to sum. m[i,j]X^(i-1)Y^(j-1);

Error, bad index into Matrix

>	And this error is due to the very same thing, since the error happens before sum receives any argument. sum(sum(m[i,j]X^(i-1)Y^(j-1),j=1..2),i=1..4);

Error, bad index into Matrix

>	The add command has special evaluation rules, which delays the evaluation of its first argument until i and j get replaced by actual numbers. add(add(m[i,j]X^(i-1)Y^(j-1),j=1..2),i=1..4);

Download Waters_01.mw

int...

Answered: acer 32405

November 08 2023

1 0

Did you include multiplication between adjacent sets of brackets?

Download ellint01.mw

note: If you use the (default) of 2-D Input mode then you can also put a space between the bracket pairs, to denote multiplication implicitly. But without either a space or a multiplication symbol it gets parsed like a function call rather than as a product of bracketed terms.

solve...

Answered: acer 32405

November 07 2023

0 6

Your phrasing is not entirely clear.

If by "...in terms of DC and T (and tau)" you mean that you you want some remaining variables expressed in terms of DC, T, and tau, then try the call,

solve({eq || (1 .. 4)}, {toff, ton, x1, x2})

If you meant something else -- eg. solving for DC,T,tau -- then I think that the phrasing was off.

Trace...

Answered: acer 32405

November 06 2023

1 0

You only need to compute a Matrix Trace once, symbolically.

Once simplified, the resulting formula can be computed pretty quickly at your sampling points.

Also, it's more efficient to use seq than it is to repeatedly augment a list.

>

restart;

>	with(LinearAlgebra):

>

>	f := x -> exp(2IPix): S1 := (x, y) -> Matrix(3, 3, [[f(x - y), 0, 0], [0, 1, 0], [0, 0, f(-y)]]): S2 := (x, y) -> Matrix(3, 3, [[f(y), 0, 0], [0, f(y - x), 0], [0, 0, 1]]): C := 1/sqrt(3)Matrix(3, 3, [[1, 1, 1], [1, f(1/3), f(2/3)], [1, f(2/3), f(1/3)]]):

>	d := Determinant(C . C);

>	V := (x, y) -> (((S1(x, y)/d^(1/3)) . C) . (S2(x, y))) . C:

>	T:=unapply(simplify(evalc(Trace(V(x,y)))),x,y);

proc (x, y) options operator, arrow; (1/6)*(-I*3^(1/2)-1)*cos(2*Pi*(x-y))+(1/6)*(3*I+3^(1/2))*sin(2*Pi*(x-y))+(1/6)*(-I*3^(1/2)-1)*cos(2*Pi*y)+(1/6)*(3*I+3^(1/2))*sin(2*Pi*y)-((1/6)*I)*3^(1/2)+(1/3)*3^(1/2)*sin(2*Pi*x)+(1/3)*cos(2*Pi*x)+1/2 end proc

>	str := time[real](): tracevals:=[seq(seq(evalhf(T(0.01i, 0.01j)), i=1..100), j=1..100)]: time[real]()-str;

>	plots:-complexplot(tracevals, x = -2 .. 3, y = -3 .. 3);

>	plots:-display( plottools:-transform((x,y,z)->[x,y])(plot3d([Re(T(x,y)), Im(T(x,y)), 1], x = 0 .. 1, y = 0 .. 1)), view=[-3..3,-3..3]);

>

Andreas132_1.mw

Ok, I doubled up calls to T in the second plot, but you could even make that faster through memoization.

I suppose that there are other ways to get the 2D plot.

one way...

Answered: acer 32405

November 06 2023

3 0

This is not terribly slow.

Check for correctness. This is a purely numeric approach. I didn't consider whether there is something smarter (eg. parametrization of the eigenvalues).

If the root-surfaces passed through each other (not happening here) then the coloring would still be by vertical position, rather than by surface.

>

restart;

>	with(LinearAlgebra):

>	f := x -> exp(2IPi*x):

>	S1 := (x, y) -> Matrix(3, 3, [[f(x - y), 0, 0], [0, 1, 0], [0, 0, f(-y)]]): S2 := (x, y) -> Matrix(3, 3, [[f(y), 0, 0], [0, f(y - x), 0], [0, 0, 1]]):

>	C := 1/sqrt(3)*Matrix(3, 3, [[1, 1, 1], [1, f(1/3), f(2/3)], [1, f(2/3), f(1/3)]]):

>	d := Determinant(C . C);

>	V := (x, y) -> (((S1(x, y)/d^(1/3)) . C) . (S2(x, y))) . C:

>	M := simplify(V(x,y)):

>	bands := proc(x, y) option remember; if not [x,y]::list(numeric) then return 'procname'(args); end if; sort(simplify(fnormal(-I*ln~(Eigenvalues(evalf(eval(M,[':-x'=x,':-y'=y]))))), 'zero')); end proc:;

>	plot3d({bands(x,y)[1],bands(x,y)[2],bands(x,y)[3]}, x=0..1, y=0..1, color=[red,cyan,blue]);

Download eigenplotthing01.mw

unevaluated return from call to h...

Answered: acer 32405

November 05 2023

6 6

The purpose of the special code about which your's asking is to have a call to procedure h return unevaluated (ie. as the same call, itself) in the case that both its arguments are not numeric.

Put briefly, it was done so that the call h(a,b) simply returns as itself, unevaluated, when a and b are merely unassigned, symbolic names rather than actual numeric values.

If h is not set up that way then the calls like h(a,b) -- that get invoked by the call BC2(a,b) -- would induce an error in dsolve/numeric.

Your call to implicitplot contains the call BC2(a,b) in the list that is passed as its first argument. That will get evaluated up-front -- in Maple's usual manner of handling calls to procedures/commands -- before implicitplot starts replacing the names a and b with actual numbers.

Below I show you that -- with your original variant -- the call BC2(a,b) returns with unevaluated calls like h(a,b) (ie. with a and b as symbolic names as its arguments).

Without the special code part of h (about which you're asking) that call BC2(a,b) would throw an error (because h(a,b) in turn would throw an error) before implicitplot even saw its arguments or had any opportunity to start replacing names a and b with actual numbers.

But first I'll show immediately below a variant in which that special construction of h is not needed. I get around the difficulties by using a so-called operator-form calling sequence to the implicitplot command. In this case implicitplot's argument consists of operators/procedures, and BC2 does not get called until its arguments are actually supplied by -- implicitplot's internals -- as numeric values.

>

restart;

>	ode := diff(y(x),x,x) + (y(x))^2*diff(y(x),x) + cos(y(x)) = x:

>	h:=proc(a,b) local F; Digits := 10; F:=dsolve({ode,y(0)=a,D(y)(0)=b},y(x),numeric); [rhs(F(1)[2]),rhs(F(1)[3])]; end proc:

>	BC1 := proc(a,b) (1 + a)b - 3a end proc: BC2 := proc(a,b) (1 + h(a,b)[1])h(a,b)[2] - 3h(a, b)[1] - h(a,b)[1]^2 + 0.5 end proc:

>	# This call produces an error. h(a,b);

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

Error, (in F) parameter 'a' must be assigned a numeric value before obtaining a solution

>	# # This call produces an error, because it calls h(a,b) # which we just saw as throwing an error. # # But we can avoid using that call BC2(a,b) (with symbolic # arguments `a` and `b`) by using an operator-form # calling sequence in the first argument to `implicitplot`. # BC2(a,b);

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

Error, (in F) parameter 'a' must be assigned a numeric value before obtaining a solution

>

#
# In this variant we don't pass BC2(a,b) directly
# to implicitplot. Instead we wrap that in an operator.
#
# Notice that the 2nd and 3rd arguments are now pure ranges,
# and not a=range, b=range
#

plots:-implicitplot([(a,b)->BC1(a,b),(a,b)->BC2(a,b)], -5..1, -11..7,
color=[black,red], legend=["BC1(a,b)=0","BC2(a,b)=0"], size=[300,300]);

>	# # And now your original. # restart;

>	ode := diff(y(x),x,x) + (y(x))^2*diff(y(x),x) + cos(y(x)) = x:

>	h:=proc(a,b) local F; Digits := 10; if type(a,numeric) and type(b,numeric) then F:=dsolve({ode,y(0)=a,D(y)(0)=b},y(x),numeric); [rhs(F(1)[2]),rhs(F(1)[3])]; else 'h'(a,b); end if; end proc:

>	BC1 := proc(a,b) (1 + a)b - 3a end proc: BC2 := proc(a,b) (1 + h(a,b)[1])h(a,b)[2] - 3h(a, b)[1] - h(a,b)[1]^2 + 0.5 end proc:

>	# # This call produces no error, since here a call to `h` # is defined to return unevaluated when its arguments # are not both numeric (eg. symbolic names). # # So we can use this expression directly in the first # argument passed to implicitplot. # BC2(a,b);

>	plots:-implicitplot([BC1(a,b),BC2(a,b)], a=-5..1, b=-11..7, color=[black,red], legend=["BC1(a,b)=0","BC2(a,b)=0"], size=[300,300]);

Download Mikey_01.mw

recovered file...

Answered: acer 32405

November 04 2023

3 1

The following attachment is as much as I believe can be recovered from your zipped .bak file.

STAT_131_Introduction_to_Prob_acc.mw

There might be two missing (unrecoverable) Inputs, in section "Day 14", between the line with "bivariate random vector" and the line with "Joint distribution functions".

I suggest that you upgrade to Maple 2023, if available to you. It seems to have fewer incidences of this kind of problem.

two ways...

Answered: acer 32405

November 02 2023

1 2

You can make the animation from your original "Test" collection, by passing the extra option insequence to your call to the plots:-display command.

Or you can account for the fact that plots:-animate passes a float value for i, by turning that into an integer prior to the indexing into your Matrix data.

Note. By default animate uses 25 frames. You could also supply the option frames=18 if you want animate to generate exactly 18 frames for i=1..18. That would make the two approaches produce the same animation. But it'd still pass i as a float, so you'd still need to trunc/floor/etc to get an integer. This option also provides an easy way to slow it down -- ie. duplicating frames by using a higher frames value, eg. multiples of 18 here, say.