787 Reputation

9 years, 146 days

Funny...

A few workarounds

b := add(rand(1 .. 6)(), i=1..3):
b := +(seq(rand(1 .. 6)(), k=1..3)):


Using the Statistics package makes the things clearer

 > restart:
 > with(Statistics):
 > Dice_1 := RandomVariable(EmpiricalDistribution([$1..6])); Dice_2 := RandomVariable(EmpiricalDistribution([$1..6])); Dice_3 := RandomVariable(EmpiricalDistribution([$1..6]));  (1)  > Sample(Dice_1, 10) +~ Sample(Dice_2, 10) +~ Sample(Dice_3, 10)  (2)  > # Or better DiceSum := Dice_1 + Dice_2 + Dice_3; Sample(DiceSum, 10)  (3)  > # what rand(1..6)() + rand(1..6)() + rand(1..6)() does # is basically Sample( Dice_1 + Dice_1 + Dice_1, 10);  (4)  > Download Statistical_equivalence.mw One way... eq1 := 2*A*B/(A*m+m)+C^2+D^2; SET := A=1/alpha: eval(eq1, SET): map(simplify, %): eval(%, isolate(SET, alpha)); 2 B 2 2 ----------- + C + D /1 \ m | - + 1| \A /  One_Way.mw Maybe this?... The matrix M depends on several quantities: • lambda • m0 • a • c • N One could write a procedure which build M from these five quantities and computes its characteristic polynomial wrt some auxiliary quantity eta. I didn't go in this direction and kept the things simple. In the attached file you will find: • a quite simple way to build M for arbitrary values of N, • the computation of its characteristic polynomial thanks to LinearAlgebra:-CharacteristicPolynomial, • the computation of this same characteristic polynomial by explicitely computing the determinant ofmatrix M-eta*Identity(N+1) wrt the last row of this matrix. CP.mw Confusing notation... Your notation FT{d(U(x, t)/dt} = d(FT{U(t))} / dt, is confusing (see @Carl Love for another interpretation) • if U(x, t) represents a function in the physical space, then its Fourier transform is: • some function H(xi, t) if you do a Fourier transformation wrt x, • some function H(x, tau) if you do a Fourier transformation wrt t. • So do we read the rhs of your equation H(xi, t) or H(x, tau) ? Given you write d(FT{U(t))} / dt, I'm leaning to diff(H(xi, t), t) (otherwise you should have written diff(H(x, tau),tau) ), meaning the Fourier transform is likely taken wrt x: fourier : U(x, t) --> H(xi, t) restart: with(inttrans): Notations: A function named by a lower case letter applies in the physical space. The same letter capitalized denotes the fourier fransform of a "physical" function. I_set := fourier(u(x, t), x, xi) = U(xi, t); fourier(u(x, t), x, xi) = U(xi, t) We_have := 'fourier'(diff(u(x, t), t), x, xi) = fourier(diff(u(x, t), t), x, xi); / d \ d fourier|--- u(x, t), x, xi| = --- fourier(u(x, t), x, xi) \ dt / dt By_notation := rhs(We_have) = eval(rhs(We_have), I_set); d d --- fourier(u(x, t), x, xi) = --- U(xi, t) dt dt Then := 'fourier'(diff(u(x, t), t), x, xi) = rhs(By_notation) / d \ d fourier|--- u(x, t), x, xi| = --- U(xi, t) \ dt / dt  Note that a Fourier transform wrt t gives I_set := fourier(u(x, t), t, tau) = U(x, tau); fourier(u(x, t), t, tau) = U(x, tau) 'fourier(diff(u(x, t), t), t, tau)' = eval(fourier(diff(u(x, t), t), t, tau), I_set); / d \ fourier|--- u(x, t), t, tau| = I tau U(x, tau) \ dt /  ... the rhs has nothing to do with diff(U(x, tau), tau) ! fourier.mw An idea... restart: kernelopts(version) Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895 a := 2*m: b := 3*m: # This can be done differently in more recent versions 'a*b' = (a)*(b); a b = (2 m) (3 m) (a)*(b) = a*b 2 (2 m) (3 m) = 6 m a*b := (a)*(b) = a*b 2 a*b := (2 m) (3 m) = 6 m  Solving numerically the Lane-Emden equat... Solving numerically the Lane-Emden equation with theta__0(xi) being itself the numerical approximation of the Emden equation returns this error message Error, (in dsolve/numeric/make_proc) ode system is singular at the initial point  Note that this has nothing to do with the fact the Emden solution is solved numerically: EmdenE := unapply( rhs(dsolve({ode1(1), theta__0(0) = 1, (D(theta__0))(0) = 0}, theta__0(xi))), xi); ode2E := n -> diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = -EmdenE(xi)^n*xi^2; ic2 := gamma(0) = 0, D(gamma)(0) = 0; LaneEmdenN := dsolve({ode2E(1), ic2}, numeric, parameters=[alpha]) Error, (in dsolve/numeric/make_proc) ode system is singular at the initial point  The reason why dsolving/numeric the Lane-Emden equation comes probably from the nullity at xi=0 of the denominator of its exact solution: LaneEmdenE := dsolve({ode2E(1), ic2}) (cos(xi alpha) xi alpha - sin(xi alpha)) _C2 gamma(xi) = -------------------------------------------- xi / / 2 \ 2\ \2 + \alpha - 1/ xi / sin(xi) - 2 xi cos(xi) + --------------------------------------------- 2 2 xi (alpha - 1) (alpha + 1)  See at the end of this answer how to et the mw file. So the idea to rewrite the initial conditions this way ic2 := gamma(epsilon) = 0, D(gamma)(epsilon) = 0;  where epsilon is aimed to be a small positive quantity (please to the first point in Remarks above). The mw file will be attached as soon as this will be possible (the site experiences a technical problem presently). Meanwhile here ifs the code (the important lines whowing how to use the numerical solution of the Emden equation are yellow highlighted) restart local gamma: ode1 := n -> (diff(xi^2*(diff(theta__0(xi), xi)), xi))/xi^2 = -theta__0(xi)^n: EmdenN := dsolve({ode1(1), theta__0(0) = 1, (D(theta__0))(0) = 0}, theta__0(xi), numeric): plots:-odeplot(EmdenN, [xi, theta__0(xi)], xi=0..1): f := proc(z) if not type(evalf(z),'numeric') then 'procname'(z); else eval(theta__0(xi), EmdenN(z)); end if; end proc: ode2 := n -> diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = -f(xi)^n*xi^2: ic2 := gamma(epsilon) = 0, D(gamma)(epsilon) = 0: LaneEmdenN := dsolve(eval({ode2(1), ic2}, epsilon=0), known=f, numeric, parameters=[alpha]) Error, (in dsolve/numeric/make_proc) ode system is singular at the initial point # define an extra parameter epsilon LaneEmdenN := dsolve({ode2(1), ic2}, known=f, numeric, parameters=[alpha, epsilon]): # with epsilon=0 LaneEmdenN(parameters=[2.356, 0]): plots:-odeplot(LaneEmdenN, [xi, gamma(xi)], xi=0..10): Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up # with epsilon=1e-6 LaneEmdenN(parameters=[2.356, 1e-6]): plots:-odeplot(LaneEmdenN, [xi, gamma(xi)], xi=1e-6..10): # works well  Remarks: • In your file you used 2 conditions at xi=0 and another one at xi=xiI: that's one to many. • I used ode1 := n -> ... (the same for ode2) to give you the possibility to change in a simple way the value of n. • I used the 'parameters' option for the same flexibility reasons. Of course you can set n, alpha and epsilon to numerical values BEFORE dsolving the equations: data := {n = 1, alpha = 2.356, epsilon=1e-6}: ode1 := (diff(xi^2*(diff(theta__0(xi), xi)), xi))/xi^2 = -theta__0(xi)^n: ic1 := theta__0(0) = 1, (D(theta__0))(0) = 0: EmdenN := dsolve(eval({ode1,ic1}, data), numeric); ode2 := diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = f(xi)^n*xi^2; ic2 := gamma(epsilon) = 0, D(gamma)(epsilon) = 0; f := proc(z) if not type(evalf(z),'numeric') then 'procname'(z); else eval(theta__0(xi), EmdenN(z)); end if; end proc; LaneEmdenN := dsolve(eval({ode2(1), ic2}, data), known=f, numeric)  Last but not least, You will see in the mw file that the exact and numerical solutions of the Lane-Emden solutions are quite different. I haven't been able to fix this in a satisfactory way. To overcome the present technical difficulties Mapleprimes faces, you can download the Maple file this way: • edit this link DropBox, • copy/paste it on your navigator, • download the file (LaneEmdenNumeric.txt), • change its name into LaneEmdenNumeric.mw, • open LaneEmdenNumeric.mw with Maple. Lazzy as I am I would try to use what ot... Maybe this could help For odd orders only: The site is in French but the Maple code is in English :-) http://villemin.gerard.free.fr/Wwwgvmm/CarreMag/aaaMaths/Tapis.htm For odd orders and few even orders (Still in french) Go here http://carres-magiques.com/partie-3-nouvelle Then search Annexe 23 and click on "Cliquez ici pour obtenir les programmes Maple" At last unzip the package This way... s := [{f0(t, r) = 0, n(t) = n(t)}, {f0(t, r) = 0, n(t) = 0}, {f0(t, r) = 0, n(t) = 0}]: pattern := n(t)=n(t): select(has, s, pattern) [{f0(t, r) = 0, n(t) = n(t)}] pattern := n(t)=0: select(has, s, pattern) [{f0(t, r) = 0, n(t) = 0}, {f0(t, r) = 0, n(t) = 0}]  For this kind question randpoly is a use... For this kind question randpoly is a useful tool. But you can also write your own generator to create equations of a specific form. For instance  > restart  > randomize()  (1)  > # your first example randpoly(x, degree=4, coeffs=rand(-5..5))  (2)  > # your second and third examples # deg := maximum degree of the polynom # E := max value of the exponent of the polynom # C := max value of the coefficients of the polynom # # radfavor = 1 favors radicals # radfavor = 0 unfavors radicals # other values of radfavor are neutral gen := proc(deg, E, C, radfavor) local a := rand(0..1): local d := rand(1..deg): local e := rand(1..E): local s := if(rand(0..1)()=0, -1, 1); local c := rand(1..C): local t, eq: if radfavor=1 then t := if(rand(0..3)()<=2, -1, 1); elif radfavor=0 then t := if(rand(0..3)()=0, -1, 1); else t := if(rand(0..1)()=0, -1, 1); end if: eq := 0: while not(has(eq, x)) do eq := s*c() * add(s*c()*x^(k*a()), k=0..deg)^(e()^t()) + s*c() = 0: end do: return eq end proc: for k from 1 to 10 do gen(2, 5, 6, 1); end do;  (3)  > for k from 1 to 10 do gen(2, 5, 6, 0); end do;  (4)  > Explanation of the parameter radfavor. Replacing t := if(rand(0..3)()<=2, -1, 1); by t := if(rand(0..19)()<=18, -1, 1); will generate a radical in 95% of the cases. If you write instead t := if(rand(0..19)()=0, -1, 1); a radical will be generated only 1 out of 20 times. Have you any idea of the size of Bd?... restart B := Matrix(9$2, symbol=m):
Bd:= CodeTools:-Usage( LinearAlgebra:-Determinant(B) ):
memory used=0.82GiB, alloc change=0.60GiB, cpu time=17.64s, real time=11.98s, gc time=9.07s

# approximately the number of symbols used to represent Bd

L := length(Bd)
40642561
# printing Bd on a 400 charaters per line display will require
# this number of lines

iquo(L, 400) + if(irem(L, 400)=0, 0, 1)
101607

# And thus a pocket book (50 characters per line and 60 lines

round(evalf(L / 50 / 60))
13548


So how can you even read this expression?

I would say no...

On might think using the package ImageTools, providing the image has an appropriate format and is of good quality.
But it will be quite difficult to extract the contours of the points, take their centers, extract the correct scales (did I miss comething?)

My advice is to use a digitizer and do the things manually (some can be use online, for instance https://plotdigitizer.com/app)

Computation of the determinant of B and ...

... is to compute the determinant X of some symbolic matrix A of elements xi, j, and next to replace in the expression of X each  xi, j by the corresponding expression in B (x[i, j]=B[i, j])
The corresponding lines are in brown courier font.

EDITED:
Following the (judicious) remark @acer  made, here is a simple say to get the determinant of B1 (not the one of Omeg):

detH  := Determinant(H):
detB1 := Determinant(B1):
detB  := detH * detB1:


 >
 >
 >
 >
 >
 >
 (1)
 >
 >
 (2)
 >
 >
 >

 > t0 := time(): detH  := Determinant(H): detB1 := Determinant(B1): detB  := detH * detB1: time()-t0;
 (3)
 > simplify(detH, size); simplify(detB1, size)
 (4)
 >
 >

As matrix Omeg := B + idn8 = H.B1 + idn8 (BTW idn8 := IdentityMatrix(8)) and because H:=Matrix(8, 8, symbol=h), (thus H is formally equivalent to the matrix X used in the "trick"), I think the "trick" can be applied to compute Determinant(Omeg).
Of course you must be aware of  @acer's remark if you have to evaluate Determinant(Omeg) for some instance of H.

What loop are you talking about? A_...

UPDATED

What loop are you talking about? (shouldn't have been here)

Percentage of variation of peak heights

A_ref := eval(A, lambda = 1.3015):
x_ref := fsolve(diff(A_ref, x), x=-2..0);
y_ref := evalf(eval(A_ref, x=x_ref));

-0.7765307466
1.920953647

data_1 := [beta=0.1, Q=1.3015, lambda=0.9986];
B_1    := eval(B, data_1):
x_1    := fsolve(diff(B_1, x), x=-2..0);
y_1    := evalf(eval(B_1, x=x_1));

PercentageOfVariation := (y_1-y_ref)/y_ref*100

[beta = 0.1, Q = 1.3015, lambda = 0.9986]
-0.8959052730
1.752467429
-8.770967392

data_2 := [beta=0.3, Q=1.3015, lambda=0.9986];
B_2    := eval(B, data_2):
dB_2   := diff(B_2, x):
x_2    := fsolve(diff(B_2, x), x=-2..0);
y_2    := evalf(eval(B_2, x=x_2));

PercentageOfVariation := (y_2-y_ref)/y_ref*100

[beta = 0.3, Q = 1.3015, lambda = 0.9986]
-1.088722291
1.547856941
-19.42247313


The main problem with Word is that the c...

The main problem with Word is that the character sizes will automatically decrease as the "depth" of the continued fraction increases (I don't know if it's possible to overcome this?) :

Why not use the following notation instead (see https://en.wikipedia.org/wiki/Continued_fraction, section Motivation and notation) and just copy-paste the result into your MS-Word document (next convert it as an equation if you prefer)?

phi := convert((1+sqrt(5))*(1/2), confrac)
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
# or
phi := [phi[], #mo("&#x2026;")]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...]


See here https://www.wikihow.life/Start-Working-with-Continued-Fractions for some advices, in particular topic (4) "to avoid the cumbersome staircase notation"

PS: note that the conventional notation should be [1 ; 1, 1, ...] instead of [1 , 1, 1, ...].

https://www.mapleprimes.com/questions/236897-Differentiate-Three-Regions-With-Implicitplot#comment296558

 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 > WhyNot := proc (alpha, delta) if not [alpha, delta]::(list(numeric)) then return ('procname')(args) end if; FirmModelHmax(alpha, .2, delta)[3] end proc: pltHmax3 := plot(   [seq(WhyNot(alpha, delta),delta=0.1..0.5,0.2)]     , alpha=0..1       , linestyle=[dot,dashdot,dash]       , legend=[seq('delta'=delta,delta=0.1..0.5,0.2)]       , legendstyle=[location=left]       , labels=["alpha","SiS firm profit"]       , labeldirections =["horizontal", "vertical"]     , legendstyle=[location=bottom] ): display(pltHmax3);
 > display([pltPP3, pltFC3, pltHmax3])
 >