Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 163 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Use intersectplot()...

tomleslie 13876
June 04 2022
3 1

as in the attached.

For some reason the intersection curve does not render well (at all?) on this site, but I can assure you that it does in the worksheet

> 

  restart;
  with(plots):
  m := 1:
  n := 2:
  Zc := (x, y) -> m*cos(1/180*x*Pi) + n*cos(1/180*y*Pi);
  Zs := (x, y) -> m*sin(1/180*x*Pi) + n*sin(1/180*y*Pi);
  display( [ plot3d
             ( [ Zc(x, y), Zs(x, y) ],
               x = 0 .. 180,
               y = 0 .. 180,
               color = [red, blue],
               transparency = 0.8,
               style = surface
             ),
             intersectplot
             ( Zc(x, y),
               Zs(x, y),
               x = 0 .. 180,
               y = 0 .. 180,
               thickness = 4,
               color = black
             )
           ]
         );

proc (x, y) options operator, arrow; m*cos((1/180)*x*Pi)+n*cos((1/180)*y*Pi) end proc

  

proc (x, y) options operator, arrow; m*sin((1/180)*x*Pi)+n*sin((1/180)*y*Pi) end proc

  

 
> 

 

Download iPlot.mw

 

In your first equation...

tomleslie 13876
June 02 2022
2 3

you have the funtion p() with argument Delta*delta - sol "collecting" on this term doesn't make sense.

As a general rule if you want to collect on say products of quantities, use a "nested" collect() command as shown in the attached

> 

  restart;
#
# If I had this expression
#
  Psi[q0]*Delta*delta - Psi[d0]*p*Delta*delta/omega[0];
#
# Then I can use a "nested" collect, like this
#
  collect(collect(%, Delta), delta);
#
# However OP has this
#
  Psi[q0]*Delta*delta - Psi[d0]*p(Delta*delta)/omega[0];
#
# where the second term contains the function p(), with
# argument Delta*delta, so one cannot "collect" on this.
#
  collect(collect(%, Delta), delta);

Psi[q0]*Delta*delta-Psi[d0]*p*Delta*delta/omega[0]

  

(Psi[q0]-Psi[d0]*p/omega[0])*Delta*delta

  

Psi[q0]*Delta*delta-Psi[d0]*p(Delta*delta)/omega[0]

  

Psi[q0]*Delta*delta-Psi[d0]*p(Delta*delta)/omega[0]

 (1)
> 

 

Download coll.mw

You have a couple of choices...

tomleslie 13876
June 02 2022
4 0

essentially dpending on whether tou want to work with T0 or omega0=2*Pi/T0. See the attached.

> 

  restart;
  with(inttrans):
  delta__t:=Sum(Dirac(t-k*T__0),k=-infinity..infinity) assuming T__0>0;
  convert(delta__t, exp);
  fourier(%, t, omega) assuming T__0>0;
  convert(%, exp);
  invfourier(%, omega, t) assuming T__0>0;

Sum(Dirac(T__0*k-t), k = -infinity .. infinity)

  

(Sum(exp((2*I)*Pi*k*t/T__0), k = -infinity .. infinity))/T__0

  

2*Pi*(Sum(Dirac(-2*Pi*k+T__0*omega), k = -infinity .. infinity))

  

-(Sum(exp(I*k*T__0*omega), k = -infinity .. infinity))

  

-(Sum(Dirac(T__0*k+t), k = -infinity .. infinity))

 (1)
> 

  restart;
  with(inttrans):
  delta__t:=Sum(Dirac(t-k*2*Pi/omega__0),k=-infinity..infinity) assuming omega__0>0;
  convert(delta__t, exp);
  fourier(%, t, omega) assuming omega__0>0;
  convert(%, exp);
  invfourier(%, omega, t) assuming omega__0>0;

Sum(Dirac(-t+2*k*Pi/omega__0), k = -infinity .. infinity)

  

(1/2)*(Sum(exp(I*k*t*omega__0), k = -infinity .. infinity))*omega__0/Pi

  

omega__0*(Sum(Dirac(k*omega__0-omega), k = -infinity .. infinity))

  

Sum(exp((2*I)*Pi*k*omega/omega__0), k = -infinity .. infinity)

  

omega__0*(Sum(Dirac(2*Pi*k+t*omega__0), k = -infinity .. infinity))

 (2)
> 

 

Download DC.mw

Depending...

tomleslie 13876
June 01 2022
1 2

on what you plan on doing with this PDE, it *may* be a better idea to define the function u(x,t) in terms of real and imaginary parts, say uR(x,t)+I*uI(x,t), in which both uR(x,t) and uI(x,t) are real functions, then convert your PDE to a system of two PDEs, each of whihc contains only real functions - something like the attached.

> 

  restart;
  interface(showassumed=0):
  assume( u__R(x,t)::real, u__I(x,t)::real);
  u(x,t):=u__R(x,t)+I*u__I(x,t):
  PDE:= evalc(diff(u(x,t),t)+diff(conjugate(u(x,t)),x)):
  PDEsys:= [ coeff(PDE, I, 0),
             coeff(PDE, I, 1)
           ];

[diff(u__R(x, t), t)+diff(u__R(x, t), x), diff(u__I(x, t), t)-(diff(u__I(x, t), x))]

 (1)
> 

 

Download RandI.mw

Seems to work...

tomleslie 13876
May 30 2022
1 1

You appear to importing data automatically (somehow) - and no-one here has access to that data, so you worksheet cannot be executed.

However if I do a cut+paste job on *some* of your data, the dataframe filtering operation *seems* to work - see the attached.

I can think of reasons why it might fail for your complete example (eg a weird datatype in column 1 maybe?) - but I think the first step is for you to download/execute the attached and see what happens


 

> 

  restart;
  interface(version);
  RawData:=Matrix([["BFH",1,"H","FY",124,124,22,2.2,0,0],
                   ["H0,075",1,"H","D1",130,153,15,3,0.05,0.06],
                   ["H0,15",1,"H","D2",136,170,13.6,2.9,0.08,0.08],
                   ["H0,2",1,"H","D3",136,181,16,2.3,0.16,0.16],
                   ["BFP",1,"P","FY",124,109,18,"",0,0],
                   ["P0,15",1,"P","D1",124,119,25,"",0,0],
                   ["P0,075",1,"P","D2",124,121,22,"",0,0],
                   ["P0",1,"P","D3",118,123,25,"",0,0],
                   ["Pneg",1,"P","D4",118,147,"","",0,0],
                   ["BFC",1,"C","FY",130,130,24,"",0,0]
                ]);

`Standard Worksheet Interface, Maple 2019.2, Windows 7, November 26 2019 Build ID 1435526`

  

Matrix(%id = 18446744074370552766)

 (1)
> 

  Data := convert(RawData, DataFrame);

DataFrame(Matrix(10, 10, {(1, 1) = "BFH", (1, 2) = 1, (1, 3) = "H", (1, 4) = "FY", (1, 5) = 124, (1, 6) = 124, (1, 7) = 22, (1, 8) = 2.2, (1, 9) = 0, (1, 10) = 0, (2, 1) = "H0,075", (2, 2) = 1, (2, 3) = "H", (2, 4) = "D1", (2, 5) = 130, (2, 6) = 153, (2, 7) = 15, (2, 8) = 3, (2, 9) = 0.5e-1, (2, 10) = 0.6e-1, (3, 1) = "H0,15", (3, 2) = 1, (3, 3) = "H", (3, 4) = "D2", (3, 5) = 136, (3, 6) = 170, (3, 7) = 13.6, (3, 8) = 2.9, (3, 9) = 0.8e-1, (3, 10) = 0.8e-1, (4, 1) = "H0,2", (4, 2) = 1, (4, 3) = "H", (4, 4) = "D3", (4, 5) = 136, (4, 6) = 181, (4, 7) = 16, (4, 8) = 2.3, (4, 9) = .16, (4, 10) = .16, (5, 1) = "BFP", (5, 2) = 1, (5, 3) = "P", (5, 4) = "FY", (5, 5) = 124, (5, 6) = 109, (5, 7) = 18, (5, 8) = "", (5, 9) = 0, (5, 10) = 0, (6, 1) = "P0,15", (6, 2) = 1, (6, 3) = "P", (6, 4) = "D1", (6, 5) = 124, (6, 6) = 119, (6, 7) = 25, (6, 8) = "", (6, 9) = 0, (6, 10) = 0, (7, 1) = "P0,075", (7, 2) = 1, (7, 3) = "P", (7, 4) = "D2", (7, 5) = 124, (7, 6) = 121, (7, 7) = 22, (7, 8) = "", (7, 9) = 0, (7, 10) = 0, (8, 1) = "P0", (8, 2) = 1, (8, 3) = "P", (8, 4) = "D3", (8, 5) = 118, (8, 6) = 123, (8, 7) = 25, (8, 8) = "", (8, 9) = 0, (8, 10) = 0, (9, 1) = "Pneg", (9, 2) = 1, (9, 3) = "P", (9, 4) = "D4", (9, 5) = 118, (9, 6) = 147, (9, 7) = "", (9, 8) = "", (9, 9) = 0, (9, 10) = 0, (10, 1) = "BFC", (10, 2) = 1, (10, 3) = "C", (10, 4) = "FY", (10, 5) = 130, (10, 6) = 130, (10, 7) = 24, (10, 8) = "", (10, 9) = 0, (10, 10) = 0}), rows = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], columns = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

 (2)
> 

  Data[(Data[1]=~"BFH")];

DataFrame(Vector[row](10, {(1) = "BFH", (2) = 1, (3) = "H", (4) = "FY", (5) = 124, (6) = 124, (7) = 22, (8) = 2.2, (9) = 0, (10) = 0}), rows = [1], columns = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])

 (3)
> 

 


 

Download DFTest.mw

Boundary conditions don't make sense...

tomleslie 13876
May 29 2022
0 7

You have the boundary conditions

 con[1][0]:=f[0](0)=(h(x)/64),(D(f[0]))(0)=0:
 con[2][0]:=g[0](0)=-(k-h(x)^2/4),(D(g[0]))(0)=0:
 for j from 1 to L do:
 con[1][j]:=f[j](0)=h(x),(D(f[j]))(0)=0:
 con[2][j]:=g[j](0)=h(x),(D(g[j]))(0)=0:
end do;

f[0](0) is the value of the function f[0](), evaluated at the point when its argument is zero. h(x) is a function defined for all 'x'. Thus the statement

f[0](0)=(h(x)/64);

is meaningless. Maybe you meant

f[0](0)=(h(0)/64);  ????

Although if this is what you mean, then it might be useful if the function h(x) were defined somehere so that h(0) can be calculated, otherwise h(0) will be treated as an in th e same way as any other simple symbol

A similar problem occurs for all other boundary conditions, sine they all contain h(x)

Fixed some typos...

tomleslie 13876
May 29 2022
2 15

in the attached.

A warning - to save execution time I only did the four  "corner" points of the Array - even this took noticeable time ~15 secs on my machine. Your original code suggests that you want to perform this calculation approximately 1million times.

You will be waiting for a loonnngggg time!

> 

  restart:
  doCalc:= proc( xi , u)  #u is the \bar(H): normalize magnetic field magnitude,
                          # where H = bar(H)*H__a
                 # Import Packages
                 uses ArrayTools, Student:-Calculus1, LinearAlgebra,
                      ListTools, RootFinding, plots, ListTools:
                 local gamma__1:= .1093,
                       alpha__3:= -0.1104e-2,
                       k__1:= 6*10^(-12),
                       d:= 0.2e-3,
                       theta0:= 0.1e-3,
                       eta__1:= 0.240e-1, chi:= 1.219*10^(-6),
                       alpha:= 1-alpha__3^2/(gamma__1*eta__1),
                       theta_init:= theta0*sin(Pi*z/d),
                       H__a:= Pi*sqrt(k__1/chi)/d,
                       H := u*H__a,     
                       c:=alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2)),
                       w := chi*H^2*eta__1*alpha/(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2),
                       n:= 20,
                       g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
                       qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
                       soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
                       AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes

                 g:= q-(1-alpha)*tan(q)- (w*q + c)*tan(q):
                 f:= subs(q = x+I*y, g):
                 b1:= evalc(Re(f)) = 0:
                 b2:= evalc(Im(f)) = 0:
                 qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
                 OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);

# Compute Odd asymptote

                 ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));
                 qstarTemporary:= min(ModifiedOddAsym);
                 indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
                 qstar2:= OddAsymptotes(indexOfqstar2);

# Compute Odd asymptote

                 AreThereComplexRoots:= type(true, 'truefalse');
                 try
                      soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});
                      soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});
                      qcomplex1:= subs(soln1, x+I*y);
                      qcomplex2:= subs(soln2, x+I*y);
                 catch:
                       AreThereComplexRoots:= type(FAIL, 'truefalse');
                 end try;

# Compute the rest of the Roots

                 OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
                 AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));
                 if   AreThereComplexRoots
                 then gg:= [ qcomplex1,
                             qcomplex2,
                             op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
                             seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)
                          ];
                 elif not AreThereComplexRoots
                 then gg:= [ op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
                             seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)
                           ];
                 end if:

# Remove the repeated roots if any & Redefine n

                 qq:= MakeUnique(gg):
                 m:= numelems(qq):

## Compute the `τ_n`time constants

                 for i to m do
                 p[i] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2 - alpha__3*xi/eta__1- chi*H^2):
                 end do:

## Compute θ_n from initial conditions

                for i to m do
                Efun[i] := cos(qq[i]-2*qq[i]*z/d)-cos(qq[i]):
                end do:

## Compute integral coefficients for off-diagonal elements θ_n matrix

                printlevel := 2:
                for i to m do
                    for j from i+1 to m do
                        b[i, j] := int(Efun[i]*Efun[j], z = 0 .. d):
                        b[j, i] := b[i, j]:
                        aa[i, j] := b[i, j]:
                        aa[j, i] := b[j, i]:
                    end do :
                end do:

## Calculate integral coefficients for diagonal elements in theta_n matrix

                for i to m do
                aa[i, i] := int(Efun[i]^2, z = 0 .. d):
                end do:

## Calculate integrals for RHS vector

               for i to m do
               F[i] := int(theta_init*Efun[i], z = 0 .. d):
               end do:

## Create matrix A and RHS vector B

               A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]):
               B := Vector([seq(F[i], i = 1 .. m)]):

## Calculate inverse of A and solve A*x=B

              Ainv := 1/A:
              r := MatrixVectorMultiply(Ainv, B):

## Define Theta(z,t)
             theta_sol := add(r[i]*Efun[i]*exp(-t/p[i]), i = 1 .. m):

## Compute v_n for n times constant

             for i to m do
             v[i] := (-2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]) + alpha__3*chi*H^2/(2*eta__1*qq[i]*gamma__1*alpha):
             end do:

## Compute v(z,t) from initial conditions
             for i to m do
             Vfun[i] := d*sin(qq[i]-2*qq[i]*z/d)+(2*z-d)*sin(qq[i]):
             end do:

## Define v(z,t)
             v_sol := add(v[i]*r[i]*Vfun[i]*exp(-t/p[i]), i = 1 .. m):

##
             minp:=min( seq( Re(p[i]), i=1..m) ):
             member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'):

## Return all the plots
                 return theta_init, theta_sol, v_sol, minp, eval(p), nstar, p[nstar], g;

                 end proc:
> 

#Convert and store xi and doCalc into an array
#
# For the sake of speed, just do the end-points of
# the requested sequences
#
#  ans1:= Array([seq(seq( [j, doCalc(j,u)], j=-2..0, 0.006), u =1..334)]):

  ans1:= CodeTools:-Usage(Array([seq(seq( [j, doCalc(j,u)], j=-2..0, 2), u =1..334,333)])):
 

memory used=1.48GiB, alloc change=112.00MiB, cpu time=15.32s, real time=14.79s, gc time=1.15s

  
> 

p_min:= convert(ans1[..,5],list):
  plot( convert(ans1[..,1],list), p_min/~30,style=line,
        labeldirections = ["horizontal", "vertical"],
        axes=boxed,
        labels =[xi, typeset(min(Re(tau[n])))],
        labelfont = ["TimesNewRoman", 14],
        axesfont = [14, 14]
     );

 
> 

 

NULL


 

Download sProc.mw

The "solution" you provide...

tomleslie 13876
May 29 2022
3 2

does not seem to be valid for your ODE system - see the attached

> 

  restart;
> 

#
# Define system and solve
#
  sys:=[ diff(y(x),x) = y(x)^4 * cos(x) + y(x)*tan(x),
         y(0) = 0.5
       ];
  mapleSol:=simplify(dsolve( sys ));

[diff(y(x), x) = y(x)^4*cos(x)+y(x)*tan(x), y(0) = .5]

  

y(x) = (cos(x)^4*(73*cos(x)^2-9)^2*(3*sin(x)+8*cos(x)))^(1/3)/(cos(x)^2*(73*cos(x)^2-9))

 (1)
> 

#
# OP *claims* solution should be
#
  OPsol:= y(x)=1/root(cos(x)^2 * 8*cos(x)-3 * sin(x),3);
#
# Plot both solutions as a quick check
#
  plot( [rhs(mapleSol), rhs(OPsol)], x=0..1, color=[red, blue]);
#
# Use odetest() to doublecheck which solution is valid
#
  odetest(sol, sys);;
  odetest(OPsol, sys);

y(x) = 1/(8*cos(x)^3-3*sin(x))^(1/3)

  

 

[0, 0]

  

[3*sin(x)^2/(cos(x)*(8*cos(x)^3-3*sin(x))^(4/3)), 0]

 (2)
> 

 

> 

 

Download bern.mw

If you had tried reading the help...

tomleslie 13876
May 27 2022
1 1

at ?strings you would find the following

Two successive double quotes that appear after the opening of a string are parsed as a single double quote.  Thus,  "abc""de" yields the string abc"de. Alternatively, a double quote can be written within a string by preceding it with a backslash, as in "abc\"de".

Is this difficult to understand?

 

Use Units:-Split()...

tomleslie 13876
May 26 2022
2 2

as in the attached.
This particular command is very useful when dealing with commands which do not handle units

> 

  restart:
  with(Units[Simple]):
  with(plots):

  a:=15*Unit('m'):
  b:=10*Unit('m'):
  displayPoints:=pointplot([a,b]);
  displayText := textplot( [ Units:-Split(a, output = coefficient),
                             Units:-Split(b, output = coefficient),
                             "text"
                           ]
                         );

 

 
> 

 


 

Download Tplot2.mw

A workaround - maybe...

tomleslie 13876
May 26 2022
3 0

There are two attached woksheets with identical code - one of them was run in Maple 2021, the other in Maple 2022. It is fairly obvious that certain commands in the GroupTheory() package have changed between these two releases. Something as simple as what is returrned a the constructor command - sometimes it is a description and sometimes it is a list of generators.

A possible workaround for your issue in Maple 2021 is to "identify" the group in the database to which your DirectProduct() output is isomorphic, and apply the DrawSubGroupLattice() command to the former - this *seems* to work

> 

  restart;
  interface(version);
  with(GroupTheory):
  G := DirectProduct(QuaternionGroup(), CyclicGroup(3));
  DrawSubgroupLattice(G);

`Standard Worksheet Interface, Maple 2022.0, Windows 7, March 8 2022 Build ID 1599809`

  

_m732551872

  

 
> 

  id:=IdentifySmallGroup(G);
  G2:=SmallGroup(id);
  DrawSubgroupLattice(G2);

24, 11

  

_m755210368

  

 
> 

 

Download GT2022.mw

 

> 

  restart;
  interface(version);
  with(GroupTheory):
  G := DirectProduct(QuaternionGroup(), CyclicGroup(3));
  DrawSubgroupLattice(G);

`Standard Worksheet Interface, Maple 2021.2, Windows 7, November 23 2021 Build ID 1576349`

  

_m656624416

  

Error, (in .) `undefined[1]` does not evaluate to a module

  
> 

  id:=IdentifySmallGroup(G);
  G2:=SmallGroup(id);
  DrawSubgroupLattice(G2);

24, 11

  

_m859921600

  

 
> 

 

Download GT2021.mw

Version issue?...

tomleslie 13876
May 26 2022
1 1

The attached using Maple 2022 seems to work just fine. Please specify the Maple version you are using and upload a worksheet (big greeen up-arrow in the Mapleprimes toolbar) illustrating your problem

> 

  restart;
  with(GroupTheory):
  G := DirectProduct(QuaternionGroup(), CyclicGroup(3)):
  DrawSubgroupLattice(G)

 
> 

  IsNormal(CyclicGroup(2), DirectProduct(CyclicGroup(2), CyclicGroup(2)));

true

 (1)
> 

 

Download GT.mw

Some obervations...

tomleslie 13876
May 25 2022
2 1

I agree with pretty much everything Carl has said on this subject.

It is particularly importatnt to determine what happens at the "final" stage. As defined the system is "conservative": there is no dissipation, so the energy injected into the system by the iniitial impulse will "rattle around"  forever,

So knowing this, what would one expect to happen? If I apply some kind of "physical intuition", then the original impulse wouldprobably propagate forwards (and maybe backwards) through the system, untl this "original impulse" became "smoothed out" and pretty much everything in the system just kept bouncing around at the same (probably low-level), forever.

There is nothing particularly difficult about setting up and solving this system for (say) 50 stages see the attached. In the attached (I think) I have set up the system so that the right-hand end is a "brick wall" - in other words the right-hand end of the right-hand springs never move - for an N-stage system one just has to specify that x[N+1](t)=y[[N+1](t)=0, whch I have done in the attached.

This still does not expalin the fundamental question of the plots you are trying to produce. The "oscillator index" is an integer function, so I do not understand how the graphs you show have a "continuous" x-axis. What the hell is an "oscillator index" of (say) 15.333333? An "oscillator index" of 15, I understand: an "oscillator index" of 16, I understand: but an "!oscillator index" of (say) 15.333333 - I have no idea what this means!

In the attached, I have used 50 stages, and plotted a representain of the" kinetic energy" in each stage (mainly to show that this can be done). I have not attempted to address the potential energy stored in the springs at each stage. Of the three springs in each stage, two of them couple between stages, so the potential energy should probably be distributed between the two coupled stages.

Anyhow - for what it is worth, and with all the provisos given above, the attached just shows that it is perfectly possible to solve 100 coupled differential equations in Maple an plot some functions of the results. I'm obviously plotting the "wrong" function - but that is because I cannot work out what it is you want to plot :-(

> 

  restart:
  with(plots):
  N:= 50: # number of stage in the system
  params:= [ m=0.022, k__g=1467.27, k__c=2.48e9, k__e=92.21, v=0.06]:
  odesys:= [ m*diff(x[1](t), t$2)+k__g*x[1](t)=-k__c*(x[1](t)-x[2](t))^3+k__e(y[1](t)-x[1](t)),
             m*diff(y[1](t), t$2)+k__g*y[1](t)=-k__c*(y[1](t)-y[2](t))^3+k__e(x[1](t)-y[1](t)),
             seq
             ( [ m*diff(x[j](t), t$2)+k__g*x[j](t)=k__c*((x[j-1](t)-x[j](t))^3-(x[j](t)-x[j+1](t))^3)+k__e(y[j](t)-x[j](t)),
                 m*diff(y[j](t), t$2)+k__g*y[j](t)=k__c*((y[j-1](t)-y[j](t))^3-(y[j](t)-y[j+1](t))^3)+k__e(x[j](t)-y[j](t))
               ][],
               j=2..N
             ),
             x[N+1](t)=0, # what happens at the right-hand end of the system??????
             y[N+1](t)=0
           ]:
  conds:= [        x[1](0)=0, D(x[1])(0)=v, y[1](0)=0, D(y[1])(0)=0,
            seq( [ x[j](0)=0, D(x[j])(0)=0, y[j](0)=0, D(y[j])(0)=0][],
                 j=2..N
               )
          ]:
  sol:= dsolve( eval( [odesys[], conds[]], params), numeric, maxfun=0):
> 

#
# Plot *some sort* of representation of kinetic energy at each stage
#
  odeplot( sol,
           [ seq( [t, eval( (m/2)*(diff(x[j](t),t)^2+diff(y[j](t),t)^2), params)],
                  j=1..N
                )
           ],
           t=0..1.5
         );

 
> 

 

> 

 

> 

 

Download bigSys.mw

Suggestion...

tomleslie 13876
May 24 2022
1 1

What do you have set under

Tools->Options->Display->Font anti-aliasing

mine is set to Enabled. (If you change this I think you have to use "Apply Globally", then shut down Maple and restart for any changes to take effect

Fixing numerous typos...

tomleslie 13876
May 22 2022
2 0

(sometimes using guesswork to figure out your intention) I can come up with the attached - which at least computes a solution for  F[3],F[4],F[5],........

However the equation to be solved at each stage contains the undefined function M(k+1). Your original code gives no indication of what this function might be, and I can't guess.

In future please use the big green up-arrow in the Mapleprimes toolbar to upload your worksheet. It save resdponders here from a lot of tedious, error-prone, retyping.

> 

  restart;
  F[0]:=0: F[1]:=0: F[2]:=A:
  for k from -1 to 10 do
      delta[k]:=0:
  od:
> 

  for k from 0 to 10 do
      F[k+3]:= solve
               ( F[k+3]+1/(k+1)/(k+2)/(k+3)*(1/2*add(cos(beta)*delta[k-m-1]*(m+1)*(m+2)*F[m+2], m=0..k)
                 -
                 1/2*sin(beta)*add(F[k-m]*(m+1)*(m+2)*F[m+2], m=0..k)+M(k+1)*F[k+1]),
                 F[k+3]
               );
od;
       

0

  

-(1/24)*M(2)*A

  

(1/60)*sin(beta)*A^2

  

(1/2880)*M(4)*M(2)*A

  

-(1/720)*sin(beta)*M(2)*A^2-(1/12600)*M(5)*sin(beta)*A^2

  

(11/20160)*sin(beta)^2*A^3-(1/967680)*M(6)*M(4)*M(2)*A

  

(1/90720)*sin(beta)*M(4)*M(2)*A^2+(1/48384)*sin(beta)*M(2)^2*A^2+(1/362880)*M(7)*sin(beta)*M(2)*A^2+(1/6350400)*M(7)*M(5)*sin(beta)*A^2

  

-(1/17280)*sin(beta)^2*A^3*M(2)-(11/4536000)*M(5)*sin(beta)^2*A^3-(11/14515200)*M(8)*sin(beta)^2*A^3+(1/696729600)*M(8)*M(6)*M(4)*M(2)*A

  

(5/266112)*sin(beta)^3*A^4-(29/958003200)*sin(beta)*M(6)*M(4)*M(2)*A^2-(7/22809600)*sin(beta)*M(4)*M(2)^2*A^2-(1/89812800)*M(9)*sin(beta)*M(4)*M(2)*A^2-(1/47900160)*M(9)*sin(beta)*M(2)^2*A^2-(1/359251200)*M(9)*M(7)*sin(beta)*M(2)*A^2-(1/6286896000)*M(9)*M(7)*M(5)*sin(beta)*A^2

  

(401/958003200)*M(4)*M(2)*A^3*sin(beta)^2+(563/319334400)*sin(beta)^2*M(2)^2*A^3+(37/479001600)*M(7)*sin(beta)^2*M(2)*A^3+(37/8382528000)*M(7)*M(5)*sin(beta)^2*A^3+(1/14784000)*M(2)*A^3*M(5)*sin(beta)^2+(1/22809600)*M(10)*sin(beta)^2*A^3*M(2)+(1/544320000)*M(10)*M(5)*sin(beta)^2*A^3+(1/1741824000)*M(10)*M(8)*sin(beta)^2*A^3-(1/919683072000)*M(10)*M(8)*M(6)*M(4)*M(2)*A

  

-(5023/2075673600)*sin(beta)^3*A^4*M(2)-(173/1945944000)*M(5)*sin(beta)^3*A^4-(23/1132185600)*M(8)*sin(beta)^3*A^4+(23/597793996800)*sin(beta)*M(8)*M(6)*M(4)*M(2)*A^2+(17/19926466560)*sin(beta)*M(2)^2*A^2*M(6)*M(4)+(1/948879360)*sin(beta)*M(4)^2*M(2)^2*A^2-(5/456648192)*M(11)*sin(beta)^3*A^4+(29/1643933491200)*M(11)*sin(beta)*M(6)*M(4)*M(2)*A^2+(7/39141273600)*M(11)*sin(beta)*M(4)*M(2)^2*A^2+(1/154118764800)*M(11)*M(9)*sin(beta)*M(4)*M(2)*A^2+(1/82196674560)*M(11)*M(9)*sin(beta)*M(2)^2*A^2+(1/616475059200)*M(11)*M(9)*M(7)*sin(beta)*M(2)*A^2+(1/10788313536000)*M(11)*M(9)*M(7)*M(5)*sin(beta)*A^2

 (1)
> 

 

Download loopSolve.mw

First 24 25 26 27 28 29 30 Last Page 26 of 207