Ronan

190 Reputation

6 Badges

6 years, 313 days
East Grinstead, United Kingdom

MaplePrimes Activity


These are answers submitted by Ronan

Try this, I extended your list,

L1 := [{3, 5}, {4, 5}, {4, 8, 9}, {-7, 2, 3}]:
      
L2 := {};  #empty set
                            L2 := {}
for i to nops(L1) do
L2 := L1[i] union L2
end do;
                          L2 := {3, 5}
                        L2 := {3, 4, 5}
                     L2 := {3, 4, 5, 8, 9}
                  L2 := {-7, 2, 3, 4, 5, 8, 9}
 

Try   implicitplot(cos(x)*cosh(y) = 1, x = -3 .. 3, y = -5 .. 5, gridrefine=2).

You can set higher values for gridrefine   3,4,5

http://www.yorku.ca/marko/ComPhys/Euler/Euler.html

 

you would need to do your own animation though.

Also here in a question, I posted at the time @Rouben Rostamian    gave a quaternion solution.

https://www.mapleprimes.com/questions/221298-I-Am-Looking-For-A-rotate-Type-Command#answer236768

This is an animation I did from following the blog I mentioned in the previous reply. It shows the momentum vectors.


I changed your c:= to an "if then else end if".That stops producing the Error  But I could have intrepited your statement  incorrectly.

 

restart;

NULL

fdsolve := proc({gamma:=NULL, rho:=NULL, mu:=NULL, omega:=NULL, t0:=NULL,
                 t1:=NULL, x0:=NULL, y0:=NULL,
                 N:=NULL}, params)
    local t, h, c, b, x, y, L, n, l, X, Y, f, g, s1, s2, s3, s4, s5, s6, s7, s8, s9;
    eval(F(t,x,y), params);
    f := unapply(%, [t,x,y]);
    eval(G(t,x,y), params);
    g := unapply(%, [t,x,y]);
    L := floor(1/mu);
    h := (t1 - t0)/N;
s1:= sum(((omega*(h*(n-1))^rho)^(j)*GAMMA(gamma+j))/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100);
s2:= sum(((omega*(h*n)^rho)^(j)*GAMMA(gamma+j))/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100);
s3:= sum(((omega*(h*n)^rho)^(j)*GAMMA(gamma+j))/(GAMMA(j*rho+mu+1)*factorial(j)*GAMMA(gamma)), j = 0 .. 100);
s4:= sum(((omega*(h*(n-i+1))^rho)^(j)*GAMMA(gamma+j))/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100);
s5:= sum(((omega*(h*(n-i-1))^rho)^(j)*GAMMA(gamma+j))/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100);
s6:= sum(((omega*(h*(n-i))^rho)^(j)*GAMMA(gamma+j))/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100);
s7:= sum(((omega*(h*(n-i))^rho)^(j)*GAMMA(gamma+j))/(GAMMA(j*rho+mu+1)*factorial(j)*GAMMA(gamma)), j = 0 .. 100);
s8:= sum(((omega*(h*(n-i-1))^rho)^(j)*GAMMA(gamma+j))/(GAMMA(j*rho+mu+1)*factorial(j)*GAMMA(gamma)), j = 0 .. 100);
s9:=sum(((omega*h^rho)^(j)*GAMMA(gamma+j))/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100);    
    c := (i,n) ->
        if i=0 then
              (h^mu)*(  ( ((n-1)^(mu+1))*s1 )-((n^(mu+1))*s2)+((n^(mu))*s3)  )else
              (h^mu)*( (((n-i+1)^(mu+1))*s4)+(((n-i-1)^(mu+1))*s5)-(2*((n-i)^(mu+1))*s6) )
        end if;

b := (i,n) -> (h^mu)*( ((n-i)^(mu)*s7)-((n-i-1)^(mu)*s8) );

    t := Array(0..N, i-> (1-i/N)*t0 + i/N*t1, datatype=float[8]);
    x[0], y[0] := x0, y0;                                           
    for n from 0 to N-1 do
        X[0], Y[0] :=
            x[0] + add(b(i,n+1)*f(t[i],x[i],y[i]), i=0..n),
            y[0] + add(b(i,n+1)*g(t[i],x[i],y[i]), i=0..n);
        for l from 1 to L do
            X[l], Y[l] :=
                x[0] + add(c(i,n+1)*f(t[i],x[i],y[i]), i=0..n)
                     + (h^mu)*s9*f(t[n+1], X[l-1], Y[l-1]),
                y[0] + add(c(i,n+1)*g(t[i],x[i],y[i]), i=0..n)
                     + (h^mu)*s9*g(t[n+1], X[l-1], Y[l-1]);
        end do;
        x[n+1], y[n+1] := X[L], Y[L];
        #printf("y[%d]=%a\n", n+1, y[n+1]);
    end do;
    return Array(0..N, i -> [t[i], x[i], y[i]]);
end proc

proc (params, { N := NULL, gamma := NULL, mu := NULL, omega := NULL, rho := NULL, t0 := NULL, t1 := NULL, x0 := NULL, y0 := NULL }) local t, h, c, b, x, y, L, n, l, X, Y, f, g, s1, s2, s3, s4, s5, s6, s7, s8, s9; eval(F(t, x, y), params); f := unapply(%, [t, x, y]); eval(G(t, x, y), params); g := unapply(%, [t, x, y]); L := floor(1/mu); h := (t1-t0)/N; s1 := sum((omega*(h*(n-1))^rho)^j*GAMMA(gamma+j)/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100); s2 := sum((omega*(h*n)^rho)^j*GAMMA(gamma+j)/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100); s3 := sum((omega*(h*n)^rho)^j*GAMMA(gamma+j)/(GAMMA(j*rho+mu+1)*factorial(j)*GAMMA(gamma)), j = 0 .. 100); s4 := sum((omega*(h*(n-i+1))^rho)^j*GAMMA(gamma+j)/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100); s5 := sum((omega*(h*(n-i-1))^rho)^j*GAMMA(gamma+j)/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100); s6 := sum((omega*(h*(n-i))^rho)^j*GAMMA(gamma+j)/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100); s7 := sum((omega*(h*(n-i))^rho)^j*GAMMA(gamma+j)/(GAMMA(j*rho+mu+1)*factorial(j)*GAMMA(gamma)), j = 0 .. 100); s8 := sum((omega*(h*(n-i-1))^rho)^j*GAMMA(gamma+j)/(GAMMA(j*rho+mu+1)*factorial(j)*GAMMA(gamma)), j = 0 .. 100); s9 := sum((omega*h^rho)^j*GAMMA(gamma+j)/(GAMMA(j*rho+mu+2)*factorial(j)*GAMMA(gamma)), j = 0 .. 100); c := proc (i, n) options operator, arrow; if i = 0 then h^mu*((n-1)^(mu+1)*s1-n^(mu+1)*s2+n^mu*s3) else h^mu*((n-i+1)^(mu+1)*s4+(n-i-1)^(mu+1)*s5-2*(n-i)^(mu+1)*s6) end if end proc; b := proc (i, n) options operator, arrow; h^mu*((n-i)^mu*s7-(n-i-1)^mu*s8) end proc; t := Array(0 .. N, proc (i) options operator, arrow; (1-i/N)*t0+i*t1/N end proc, datatype = float[8]); x[0], y[0] := x0, y0; for n from 0 to N-1 do X[0], Y[0] := x[0]+add(b(i, n+1)*f(t[i], x[i], y[i]), i = 0 .. n), y[0]+add(b(i, n+1)*g(t[i], x[i], y[i]), i = 0 .. n); for l to L do X[l], Y[l] := x[0]+add(c(i, n+1)*f(t[i], x[i], y[i]), i = 0 .. n)+h^mu*s9*f(t[n+1], X[l-1], Y[l-1]), y[0]+add(c(i, n+1)*g(t[i], x[i], y[i]), i = 0 .. n)+h^mu*s9*g(t[n+1], X[l-1], Y[l-1]) end do; x[n+1], y[n+1] := X[L], Y[L] end do; return Array(0 .. N, proc (i) options operator, arrow; [t[i], x[i], y[i]] end proc) end proc

(1)

 

 


 

Download _Prab-1.mw

If you want to some home study by yourself this is very good. Called Maths Foundation series.
 
https://www.youtube.com/user/njwildberger/playlists?view=50&shelf_id=10&sort=dd
 

Hello,

Dont know if this will help you or not. I posted the link below in March.

There are different definitions of the elliptic integrals used by various software packages. Without doing a couple of nights restudy I can't offer much else. Hope it helps.

https://www.mapleprimes.com/questions/221540-How-Do-I-Speed-Up-The-Evaluation-Of

 

 

180*evalf(FindAngle(p1, p2))/Pi

evalf evaluates acrcos to a decimal number in radians. 180/Pi converts radians to degrees.

Thy increasing the number of diditas displayed to 5 million and see what is that is ok.

Yes. Dr R Lopez posted this on the change in Maple 2017 http://www.mapleprimes.com/maplesoftblog/208309-Typesetting-And-Maple-2017

The Determinant to the Equations = 0 so the system cant be solved. Maybe you have an error.
 

restart

``

with(LinearAlgebra)

[`&x`, Add, Adjoint, BackwardSubstitute, BandMatrix, Basis, BezoutMatrix, BidiagonalForm, BilinearForm, CARE, CharacteristicMatrix, CharacteristicPolynomial, Column, ColumnDimension, ColumnOperation, ColumnSpace, CompanionMatrix, CompressedSparseForm, ConditionNumber, ConstantMatrix, ConstantVector, Copy, CreatePermutation, CrossProduct, DARE, DeleteColumn, DeleteRow, Determinant, Diagonal, DiagonalMatrix, Dimension, Dimensions, DotProduct, EigenConditionNumbers, Eigenvalues, Eigenvectors, Equal, ForwardSubstitute, FrobeniusForm, FromCompressedSparseForm, FromSplitForm, GaussianElimination, GenerateEquations, GenerateMatrix, Generic, GetResultDataType, GetResultShape, GivensRotationMatrix, GramSchmidt, HankelMatrix, HermiteForm, HermitianTranspose, HessenbergForm, HilbertMatrix, HouseholderMatrix, IdentityMatrix, IntersectionBasis, IsDefinite, IsOrthogonal, IsSimilar, IsUnitary, JordanBlockMatrix, JordanForm, KroneckerProduct, LA_Main, LUDecomposition, LeastSquares, LinearSolve, LyapunovSolve, Map, Map2, MatrixAdd, MatrixExponential, MatrixFunction, MatrixInverse, MatrixMatrixMultiply, MatrixNorm, MatrixPower, MatrixScalarMultiply, MatrixVectorMultiply, MinimalPolynomial, Minor, Modular, Multiply, NoUserValue, Norm, Normalize, NullSpace, OuterProductMatrix, Permanent, Pivot, PopovForm, ProjectionMatrix, QRDecomposition, RandomMatrix, RandomVector, Rank, RationalCanonicalForm, ReducedRowEchelonForm, Row, RowDimension, RowOperation, RowSpace, ScalarMatrix, ScalarMultiply, ScalarVector, SchurForm, SingularValues, SmithForm, SplitForm, StronglyConnectedBlocks, SubMatrix, SubVector, SumBasis, SylvesterMatrix, SylvesterSolve, ToeplitzMatrix, Trace, Transpose, TridiagonalForm, UnitVector, VandermondeMatrix, VectorAdd, VectorAngle, VectorMatrixMultiply, VectorNorm, VectorScalarMultiply, ZeroMatrix, ZeroVector, Zip]

(1)

``

eq1 := x[1]-x[2] = 0;

x[1]-x[2] = 0

(2)

eq2 := -x[1]+2*x[2]-x[3] = 0;

-x[1]+2*x[2]-x[3] = 0

(3)

eq3 := -x[2]+2*x[3]-x[4] = 0;

-x[2]+2*x[3]-x[4] = 0

(4)

eq4 := -x[3]+x[4]-t = 0;

-x[3]+x[4]-t = 0

(5)

M, v := GenerateMatrix({eq1, eq2, eq3, eq4}, [x[1], x[2], x[3], x[4]])

Matrix(%id = 18446744074805789630), Vector[column](%id = 18446744074805789510)

(6)

LinearSolve(M, v)

Error, (in LinearAlgebra:-LinearSolve) inconsistent system

 

LinearAlgebra:-Rank(M)

4

(7)

GaussianElimination(M)

Matrix(%id = 18446744074805811910)

(8)

LinearAlgebra:-Determinant(M)

0

(9)

``


 

Download Eqns-1.mw


This Is how I would approach the issue. Could be a problem if l=1 or n.

restart

``

C := (sum(A[k]*x[i, k], k = 1 .. l-1)+sum(A[k]*x[i, k], k = l+1 .. n)+A[l]*x[i, l])^2-(2*(sum(A[k]*x[i, k], k = 1 .. l-1)+sum(A[k]*x[i, k], k = l+1 .. n)+A[l]*x[i, l]))*y[i]+y[i]^2

(sum(A[k]*x[i, k], k = 1 .. l-1)+sum(A[k]*x[i, k], k = l+1 .. n)+A[l]*x[i, l])^2-2*(sum(A[k]*x[i, k], k = 1 .. l-1)+sum(A[k]*x[i, k], k = l+1 .. n)+A[l]*x[i, l])*y[i]+y[i]^2

(1)

diff(C, A[l])

2*(sum(A[k]*x[i, k], k = 1 .. l-1)+sum(A[k]*x[i, k], k = l+1 .. n)+A[l]*x[i, l])*x[i, l]-2*x[i, l]*y[i]

(2)

``


Download diff_sum.mw

 

 

.Replace in solve({.....................},{c__1,c__3,c__4,x}) and it works

 


You haven't defined "p", you have "ep" so I assumed that to be e*p. You also have defined "del" so I changed that to 1 just to see what would happen. Could be totally wrong but it solves easily then.

``

restart

e := 5

5

(1)

f := 0.714e-1

0.714e-1

(2)

p := 0.815e-1

0.815e-1

(3)

q := .397

.397

(4)

Sol := solve({(-n+(p*q/((1+f*n)^2+(1+T-n)^2))^2)/(e*p) = 0, -T+p*(1+e*f*n)/((1+f*n)^2+(1+T-n)^2) = 0}, [T, n])

[[T = 0.3919624923e-1, n = 0.2421000982e-3], [T = -.9371446226-.8116256169*I, n = 0.6704759545e-1+.2196973222*I], [T = -2.934201781+3.533684884*I, n = -1.734513872+2.660302975*I], [T = -3.532762021+1.455618087*I, n = -2.704556354+2.261007132*I], [T = -3.532762021-1.455618087*I, n = -2.704556354-2.261007132*I], [T = -2.934201781-3.533684884*I, n = -1.734513872-2.660302975*I], [T = -.9371446226+.8116256169*I, n = 0.6704759545e-1-.2196973222*I]]

(5)

nops(Sol)

7

(6)

for i to nops(Sol) do Sol[i] end do

[T = 0.3919624923e-1, n = 0.2421000982e-3]

 

[T = -.9371446226-.8116256169*I, n = 0.6704759545e-1+.2196973222*I]

 

[T = -2.934201781+3.533684884*I, n = -1.734513872+2.660302975*I]

 

[T = -3.532762021+1.455618087*I, n = -2.704556354+2.261007132*I]

 

[T = -3.532762021-1.455618087*I, n = -2.704556354-2.261007132*I]

 

[T = -2.934201781-3.533684884*I, n = -1.734513872-2.660302975*I]

 

[T = -.9371446226+.8116256169*I, n = 0.6704759545e-1-.2196973222*I]

(7)

``

``


Download maple-1.mw

 


This might help you. It was posted as an answer to my animation question many months ago. I found it extremly helpful.

http://www.mapleprimes.com/posts/136156-Example-Of-Animation

 

 

Use Rank it is in the LinearAlgebra Package.

Rank(A) returns a value of   5

Hope this helps.

 

 

1 2 Page 1 of 2