Advanced Search

acer

32303 Reputation

29 Badges

19 years, 309 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

Contact acer

MaplePrimes Activity


These are answers submitted by acer

fnormal...

acer 32303
November 18 2023
3 0

I would prefer Maple's own fnormal (regardless of whether Chop were to act as you'd like).

But you could target Chop at the (complex) floats, if you'd prefer.

> 

restart;
> 

ode:=[diff(x(t), t) = -3*x(t) + 4*y(t), diff(y(t), t) = 5*x(t) + 9*z(t), diff(z(t), t) = y(t) + 6*z(t)];
sol:=dsolve(ode):

[diff(x(t), t) = -3*x(t)+4*y(t), diff(y(t), t) = 5*x(t)+9*z(t), diff(z(t), t) = y(t)+6*z(t)]

> 

temp := evalf[16](sol):
> 

simplify(fnormal(temp),zero);

{x(t) = .8172764111*c__1*exp(1.894304969*t)-1.150854760*c__2*exp(-6.475677505*t)+.3780227930*c__3*exp(7.581372536*t), y(t) = c__1*exp(1.894304969*t)+c__2*exp(-6.475677505*t)+c__3*exp(7.581372536*t), z(t) = -.2435641207*c__1*exp(1.894304969*t)-0.8015596745e-1*c__2*exp(-6.475677505*t)+.6323620634*c__3*exp(7.581372536*t)}

> 

subsindets(temp,complex(float),MmaTranslator:-Mma:-Chop);

{x(t) = .8172764111*c__1*exp(1.894304969*t)-1.150854760*c__2*exp(-6.475677505*t)+.3780227930*c__3*exp(7.581372536*t), y(t) = c__1*exp(1.894304969*t)+c__2*exp(-6.475677505*t)+c__3*exp(7.581372536*t), z(t) = -.2435641207*c__1*exp(1.894304969*t)-0.8015596745e-1*c__2*exp(-6.475677505*t)+.6323620634*c__3*exp(7.581372536*t)}

Download chop_fnormal.mw

another...

acer 32303
November 16 2023
2 1

If n>1 and m+3=n then m+3>1.

> 

restart;
> 

kernelopts(version);

`Maple 2023.1, X86 64 LINUX, Jul 07 2023, Build ID 1723669`

> 

rel := n -> (n-3)^(n/(n-1))*2^(n/(n-1))
            -((n-1)*2^(n/(n-1))-4*2^(1/(n-1)))*(n-3)^(1/(n-1)):
> 

simplify(factor(expand(rel(m+3)))) assuming m+3>1;

0

> 

simplify(factor(expand(rel(m+3)))) assuming m<>-2;

0

Download Prove_It_True_ac.mw

(Also shown above for m<>-2, ie. n<>1)

ps. I will submit a report against this weakness is simplify (as well as your original).

> 

expr := -2^((2*m + 5)/(m + 2)) + 4*2^(1/(m + 2)):
> 

simplify(expr) assuming m<>-2;

-2^((2*m+5)/(m+2))+4*2^(1/(m+2))

> 

simplify(factor(expand(expr))) assuming m<>-2;

0

Download Prove_It_True_acc.mw

W1:=W...

acer 32303
November 16 2023
2 0

Your procedure assigns W to W1.

That makes W1 be the very same Matrix as W, not a copy. The two names then point to the same Matrix.

Then the procedure acts inplace on the first column of W1, which means that W has changed.

What you probably want your procedure to do instead is,

   local W1 := copy(W);

LUDecomposition...

acer 32303
November 16 2023
3 1

Is this a correct import of your Matrix?

> 

restart;
> 

kernelopts(version);

`Maple 2023.1, X86 64 LINUX, Jul 07 2023, Build ID 1723669`

 

> 

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

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

 

> 

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

0

> 

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(%);

true

9339

-(-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)

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

Download det_fun.mw

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

acer 32303
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);

-.5578872507 .. .9544633133, -.7093652620 .. .7093652620

> 

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

view = [-.5578872507 .. .9544633133, -.7093652620 .. .7093652620], xtickmarks = [0., -.5000000000, .9000000000], ytickmarks = [0., -.7000000000, .7000000000]

 

> 

# 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...

acer 32303
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...

acer 32303
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...

acer 32303
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);

[[2, 3, 5, 1, 9], [2, 3, 7, 1, 5], [2, 3, 10, 1, 4], [2, 4, 3, 1, 6], [2, 4, 6, 1, 3], [2, 4, 9, 1, 5], [2, 4, 10, 1, 7], [2, 5, 4, 1, 8], [2, 5, 6, 1, 4], [2, 5, 7, 1, 3], [2, 6, 4, 1, 5], [2, 6, 5, 1, 3], [2, 6, 5, 1, 10], [2, 6, 8, 1, 3], [2, 6, 8, 1, 4], [2, 7, 10, 1, 4], [2, 8, 4, 1, 6], [2, 8, 6, 1, 9], [2, 8, 7, 1, 4], [2, 8, 9, 1, 7], [2, 8, 10, 1, 5], [2, 9, 4, 1, 8], [2, 9, 8, 1, 6], [2, 9, 10, 1, 4], [2, 10, 4, 1, 7], [2, 10, 9, 1, 5], [2, 10, 9, 1, 8], [3, 2, 4, 1, 8], [3, 2, 8, 1, 4], [3, 2, 8, 1, 6], [3, 4, 10, 1, 8], [3, 6, 9, 1, 5], [3, 10, 7, 1, 5], [4, 7, 10, 1, 8], [5, 1, 7, 2, 9], [5, 2, 9, 1, 7], [5, 6, 1, 2, 7], [5, 7, 6, 2, 4], [5, 7, 10, 2, 8], [5, 8, 4, 2, 3], [5, 8, 7, 2, 4], [5, 9, 8, 2, 4], [5, 10, 4, 2, 6], [5, 10, 6, 2, 3], [5, 10, 9, 2, 4], [5, 10, 9, 2, 7], [6, 4, 7, 2, 5], [6, 4, 9, 1, 7], [6, 8, 1, 2, 7], [6, 10, 5, 2, 3], [6, 10, 5, 2, 7], [6, 10, 8, 2, 4], [7, 8, 10, 1, 4], [7, 8, 10, 2, 5], [8, 4, 1, 2, 9], [8, 6, 2, 1, 3], [9, 6, 1, 2, 4], [10, 2, 1, 3, 6], [10, 2, 4, 3, 1], [10, 2, 8, 3, 7], [10, 5, 4, 3, 6], [10, 6, 9, 3, 2], [10, 7, 1, 3, 5], [10, 7, 8, 3, 2], [10, 8, 2, 3, 5], [10, 8, 7, 1, 4], [10, 8, 7, 2, 5], [10, 8, 7, 3, 2], [10, 8, 7, 3, 6], [10, 9, 6, 3, 2], [10, 9, 8, 3, 6]]

71

Download minhthien.mw

example...

acer 32303
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?

> 

restart; with(plots)
> 

Sc := 1; sigma := 5; E := 1; delta := 1; n := .5; Pr := .71; alphac := .1; alphat := .1; Lt := .1; Br := .1; Rd := 2; deltaB := .1
> 

OdeSys := {(diff(Theta(xi, eta), eta, eta))/Pr+(1/2)*(1-xi)*eta*(diff(Theta(xi, eta), eta))+xi*f(xi, eta)*(diff(Theta(xi, eta), eta))-xi*(1-xi)*(diff(Theta(xi, eta), xi)), diff(h(xi, eta), eta, eta)+(1/2)*(1-xi)*eta*(diff(f(xi, eta), eta))+xi*(f(xi, eta)*(diff(h(xi, eta), eta))-(diff(f(xi, eta), eta))*h(xi, eta)-2*lambda*(diff(f(xi, eta), eta)))-xi*(1-xi)*(diff(h(xi, eta), xi)), diff(f(xi, eta), eta, eta, eta)+(1/2)*(1-xi)*eta*(diff(f(xi, eta), eta, eta))+xi*(f(xi, eta)*(diff(f(xi, eta), eta, eta))-(diff(f(xi, eta), eta))*(diff(f(xi, eta), eta))+2*lambda*h(xi, eta))-xi*(1-xi)*(diff(diff(f(xi, eta), eta), xi)), (diff(Phi(xi, eta), eta, eta))/Sc+(1/2)*(1-xi)*eta*(diff(Phi(xi, eta), eta))+xi*f(xi, eta)*(diff(Phi(xi, eta), eta))-(sigma*sigma)*xi*Phi(xi, eta)*(1+n*delta*Theta(xi, eta))*exp(-E/(1+delta*Theta(xi, eta)))-xi*(1-xi)*(diff(Phi(xi, eta), xi))}; Cond := {Phi(0, eta) = 0, Phi(xi, 0) = 1, Phi(xi, 10) = 0, Theta(0, eta) = 0, Theta(xi, 0) = 1, Theta(xi, 10) = 0, f(0, eta) = 0, f(xi, 0) = 0, h(0, eta) = 0, h(xi, 0) = 0, h(xi, 10) = 0, (D[2](f))(xi, 0) = 1, (D[2](f))(xi, 10) = 0}; Cond := {Phi(0, eta) = 0, Phi(xi, 0) = 1, Phi(xi, 10) = 0, Theta(0, eta) = 0, Theta(xi, 0) = 1, Theta(xi, 10) = 0, f(0, eta) = 0, f(xi, 0) = 0, h(0, eta) = 0, h(xi, 0) = 0, h(xi, 10) = 0, (D[2](f))(xi, 0) = 1, (D[2](f))(xi, 10) = 0}
> 

colour := [red, green, blue, gold]
> 

lambdaVals := [1, 1.5, 2, 3]
> 

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+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:
> 

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:
> 

plotA := plots:-display(seq(Ans[k]:-plot(Theta(xi, eta), xi = 1, eta = 0 .. 10, color = colour[k]), k = 1 .. nops(lambdaVals)), linestyle = "solid", thickness = 2, size = [700, 400])

> 

display(seq(plot(TTheta(1, eta, k), eta = 0 .. 10, color = colour[k], thickness = 2), k = 1 .. nops(lambdaVals)), size = [700, 400])

> 

display(seq(plot(Ng(1, eta, k), eta = 0 .. .99, color = colour[k], thickness = 2, adaptive = false), k = 1 .. nops(lambdaVals)), size = [700, 400])

> 

display(seq(plot(Bj(1, eta, k), eta = 0 .. .99, color = colour[k], thickness = 2, adaptive = false), k = 1 .. nops(lambdaVals)), size = [700, 400])

NULL

Download error_in_plot_ac.mw

example...

acer 32303
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.

restart

A := 'sin(sqrt(a^2+b^2-2*a*b*sin(.123)))'

sin(sqrt(a^2+b^2-2*a*b*sin(.123)))

a := 2*Unit('m')

b := 3*Unit('m') 

``

Parse:-ConvertTo1D, "`%1` is not a module or member", InertForm

NULL

p(A)

0, "%1 is not a command in the %2 package", _Hold, Typesetting

NULL

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.

restart

A := 'sin(sqrt(a^2+b^2-2*a*b*sin(.123)))'

sin(sqrt(a^2+b^2-2*a*b*sin(.123)))

a := 2*Unit('m')

b := 3*Unit('m') 

NULL

> 

p := proc(A::uneval) local A1, B, C;
  uses InertForm, Typesetting;
  A1 := eval(A,2);
  B := eval(MakeInert(A1));
  B := subsindets(B,`*`,`%*`@op);
  C := eval(A);
  mrow(Typeset(A), mo("="),
  Typeset(A1), mo("="),
  subs("&sdot;"="&ast;",Display(B,':-inert'=false)), mo("="),
  Typeset(C));
 end proc:

NULL

p(A)

0, "%1 is not a command in the %2 package", _Hold, Typesetting

NULL

Download TS_multi_eq_uneval_2.mw

add...

acer 32303
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);

-11.518

Download Waters_01.mw

int...

acer 32303
November 08 2023
1 0

Did you include multiplication between adjacent sets of brackets?

> 

int(1/sqrt((a - t)*(t - b)*(t - c)*(t - d)), t = b .. a)
  assuming a>t, t>b, b>c, c>d;

2*EllipticK(((a-b)*(c-d)/((a-c)*(b-d)))^(1/2))/(a*b-a*d-b*c+c*d)^(1/2)

> 

convert(%,EllipticF);

2*EllipticF(1, ((a-b)*(c-d)/((a-c)*(b-d)))^(1/2))/(a*b-a*d-b*c+c*d)^(1/2)

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...

acer 32303
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...

acer 32303
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(2*I*Pi*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);

-1

> 

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.01*i, 0.01*j)), i=1..100), j=1..100)]:
time[real]()-str;

0.29e-1

> 

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...

acer 32303
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(2*I*Pi*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);

-1

> 

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

First 33 34 35 36 37 38 39 Last Page 35 of 336