Advanced Search

Paras31

130 Reputation

3 Badges

0 years, 169 days
Hellenic Open University
Mathematician

Social Networks and Content at Maplesoft.com

Contact Paras31
Teacher of Mathematics with a proven track record of working in education management. Proficient in Ease of Adaptation, Course Design, and Instructional Technology. Strong education professional with a Bachelor's degree with an emphasis in Mathematics from the University of Aegean. Currently, he is pursuing a Master's degree in Applied Mathematics at the Hellenic Open University, with a specific focus on Ordinary and Partial Differential equations. His enthusiasm lies in the application of mathematical models to real-world contexts, such as epidemiology and population growth. Aside from his passion for teaching, Athanasios enjoys football, basketball, and spending time with his dogs.

MaplePrimes Activity


These are replies submitted by Paras31

@acer you are right about that! I h...

Paras31 130
October 03 2024
0 0

@acer you are right about that! I had taken results from Mathematica as you can see in the attached document DNAsegment3.pdf.

However, the results are different

 

I am trying to apply conclusions from th...

Paras31 130
October 02 2024
0 0

I am trying to apply conclusions from the following paper on DNA potential to my research on proteins: Discrete Nonlinear Schro¨dinger Equation and Polygonal Solitons with Applications to Collapsed Proteins Specifically, I have differentiated the potential presented in the paper and embedded it into a discrete Klein-Gordon equation chain to see if it yields the desired results for my system.

This is the potential

and these are the parameters (soliton 2) I used

However, I am facing a couple of issues:

  1. Performance Issue: When I run the code, everything appears to work well initially, but when I attempt to evaluate snapshot(50) (or similar evaluations), Maple keeps calculating indefinitely without displaying any result. I am wondering if there is a way to make my code execute faster or optimize the calculation to get results without significant delays. I suspect the issue could be related to the computational complexity or the numerical methods Maple uses to handle the ODEs.

  2. Inconsistent Plot Results: When I try to plot using plots:-odeplot(sols, [t, U[67](t)], t = 0 .. 58, color = "green"), I get two different results upon repeated attempts. There seems to be some inconsistency, and I am not sure why this is happening. Could this be related to Maple's numerical solver, or perhaps something about the way sols is being evaluated? Any insights into why this inconsistency occurs and how to get a stable, consistent plot would be very helpful.

Do you have any suggestions on optimizing my Maple code, particularly for evaluating snapshot at higher steps, and fixing the plotting inconsistency would be greatly appreciated.

v1.mw

restart; with(plots); K := 99; L := 200; delta := 0.1999385999e-1
 

h := L/(K+1.0)

2.000000000

 (1)

`#mover(mi("δ",fontstyle = "normal"),mo("ˆ"))` := h*delta

0.3998771998e-1

 (2)

NULL

x := [seq(-(1/2)*L+h*i, i = 0 .. K+1)]

NULL

N := numelems(x)

101

 (3)

NULL

`#mover(mi("U"),mo("&rarr;"))` := `<,>`(seq(U[n](t), n = 1 .. N))

_rtable[36893488152041066852]

 (4)

 

 

odes := seq(ode(n), n = 1 .. N)

ode(1), ode(2), ode(3), ode(4), ode(5), ode(6), ode(7), ode(8), ode(9), ode(10), ode(11), ode(12), ode(13), ode(14), ode(15), ode(16), ode(17), ode(18), ode(19), ode(20), ode(21), ode(22), ode(23), ode(24), ode(25), ode(26), ode(27), ode(28), ode(29), ode(30), ode(31), ode(32), ode(33), ode(34), ode(35), ode(36), ode(37), ode(38), ode(39), ode(40), ode(41), ode(42), ode(43), ode(44), ode(45), ode(46), ode(47), ode(48), ode(49), ode(50), ode(51), ode(52), ode(53), ode(54), ode(55), ode(56), ode(57), ode(58), ode(59), ode(60), ode(61), ode(62), ode(63), ode(64), ode(65), ode(66), ode(67), ode(68), ode(69), ode(70), ode(71), ode(72), ode(73), ode(74), ode(75), ode(76), ode(77), ode(78), ode(79), ode(80), ode(81), ode(82), ode(83), ode(84), ode(85), ode(86), ode(87), ode(88), ode(89), ode(90), ode(91), ode(92), ode(93), ode(94), ode(95), ode(96), ode(97), ode(98), ode(99), ode(100), ode(101)

 (5)

 

W := proc (Un) options operator, arrow; -(1/4)*(b*d-e*q)^2/(b^2*(e+b*Un^2))-(1/4)*(8*b*c*m^2+q^2)*Un^2/b+c*Un^4 end proc; dW := D(W); dW(Un)

(1/2)*(b*d-e*q)^2*Un/(b*(Un^2*b+e)^2)-(1/2)*(8*b*c*m^2+q^2)*Un/b+4*c*Un^3

 (6)

#where each ode is defined by this equation:

"ode(n) := (U)[n] - `Delta__d`(n) + delta*(U)[n] +((b d-e q)^2 U[n] )/(2 b (b U[n]^(2) +e)^2)-( (8 b c m^2+q^2) U[n] )/(2 b)+4 c U[n]^(3) =0; "

proc (n) options operator, arrow, function_assign; (diff(`#mover(mi("U"),mo("&rarr;"))`, t, t))[n]-Delta__d(n)+`#mover(mi("&delta;",fontstyle = "normal"),mo("&circ;"))`*(diff(`#mover(mi("U"),mo("&rarr;"))`, t))[n]+(1/2)*(b*d-e*q)^2*`#mover(mi("U"),mo("&rarr;"))`[n]/(b*(b*`#mover(mi("U"),mo("&rarr;"))`[n]^2+e)^2)-(1/2)*(8*b*c*m^2+q^2)*`#mover(mi("U"),mo("&rarr;"))`[n]/b+4*c*`#mover(mi("U"),mo("&rarr;"))`[n]^3 = 0 end proc

 (7)

#and  the operator function is defined as:

"`Delta__d`(n) := {[[0,n=1 or n = N],[U[n+1] - 2*U[n] +U[n-1],otherwise]]; "

proc (n) options operator, arrow, function_assign; piecewise(n = 1 or n = N, 0, `#mover(mi("U"),mo("&rarr;"))`[n+1]-2*`#mover(mi("U"),mo("&rarr;"))`[n]+`#mover(mi("U"),mo("&rarr;"))`[n-1]) end proc

 (8)

 

#Define a sequence of initial conditions:

initPositions := seq(U[n](0) = U__0(n), n = 1 .. N); initVelocities := seq((D(U[n]))(0) = 0, n = 1 .. N)
 

#where:

"`U__0`(n) :={[[0.0,n = 1 or n = N],[a*sin((j*Pi*h*x[n]/L) ),otherwise]];"

proc (n) options operator, arrow, function_assign; piecewise(n = 1 or n = N, 0., a*sin(j*Pi*h*x[n]/L)) end proc

 (9)

 

 

b := -0.222159366e-3; c := 1.088046084; d := 0.1308858509e-2; e := 0.394423507e-3; q := -0.64844084e-3; m := 1.518466566

1.518466566

 (10)

 

a := (1/2)*Pi; j := 2

2

 (11)

 

``plot([seq([x[n], U__0(n)], n = 1 .. N)], style = point, axes = boxed, gridlines = true, labels = ['x__n', 'U__n(t = 0)'], labelfont = [Times, 16], size = [500, 300], color = "red")

  

sols := dsolve({odes, initPositions, initVelocities}, numeric, maxfun = 0, method = rkf45, output = listprocedure, abserr = 0.1e-8, stepsize = 0.5e-3)

`[Length of output exceeds limit of 1000000]`

 (12)

plots:-odeplot(sols, [t, U[67](t)], t = 0 .. .58, color = "green")

  

 

plots:-odeplot(sols, [t, U[67](t)], t = 0 .. 100, color = "green")

  

`#mover(mi("Usol"),mo("&rarr;"))` := `<,>`(seq(eval(U[n](t), sols), n = 1 .. N))

 

`#mover(mi("Usol"),mo("&rarr;"))`[100](3.0)

HFloat(2.985185341413548)

 (13)

" SnapshotPlot(t) := plot([seq([x[n], Usol[n](t)], n = 1..N)], style = point, symbolsize = 10, symbol = solidcircle, axes = boxed, gridlines = true, labels = ['`x__n`', '`U__n`'], labelfont = [Times, 15], title = typeset(`t = `, t), titlefont = [Garmond, 20], size = [600, 300],color="blue" ):"

Warning, (in SnapshotPlot) `n` is implicitly declared local

  

SnapshotPlot(100)

  

plots:-display(Array(1 .. 2, 1 .. 2, `<,>`(`<|>`(SnapshotPlot(10), SnapshotPlot(20)), `<|>`(SnapshotPlot(30), SnapshotPlot(3000)))))

Error, (in `#mover(mi("Usol"),mo("&rarr;"))`[1]) cannot evaluate the solution further right of 197.51949, probably a singularity

  

 

Download v1.mw

@acer I probably deleted it when I ...

Paras31 130
September 30 2024
0 0


 

restart; with(plots); K := 99; L := 200; delta := 0.5e-1
 

h := L/(K+1.0)

2.000000000

 (1)

`#mover(mi("&delta;",fontstyle = "normal"),mo("&circ;"))` := h*delta

.1000000000

 (2)

NULL

x := [seq(-(1/2)*L+h*i, i = 0 .. K+2)]

[-100., -98.00000000, -96.00000000, -94.00000000, -92.00000000, -90.00000000, -88.00000000, -86.00000000, -84.00000000, -82.00000000, -80.00000000, -78.00000000, -76.00000000, -74.00000000, -72.00000000, -70.00000000, -68.00000000, -66.00000000, -64.00000000, -62.00000000, -60.00000000, -58.00000000, -56.00000000, -54.00000000, -52.00000000, -50.00000000, -48.00000000, -46.00000000, -44.00000000, -42.00000000, -40.00000000, -38.00000000, -36.00000000, -34.00000000, -32.00000000, -30.00000000, -28.00000000, -26.00000000, -24.00000000, -22.00000000, -20.00000000, -18.00000000, -16.00000000, -14.00000000, -12.00000000, -10.00000000, -8.00000000, -6.00000000, -4.00000000, -2.00000000, 0., 2.0000000, 4.0000000, 6.0000000, 8.0000000, 10.0000000, 12.0000000, 14.0000000, 16.0000000, 18.0000000, 20.0000000, 22.0000000, 24.0000000, 26.0000000, 28.0000000, 30.0000000, 32.0000000, 34.0000000, 36.0000000, 38.0000000, 40.0000000, 42.0000000, 44.0000000, 46.0000000, 48.0000000, 50.0000000, 52.0000000, 54.0000000, 56.0000000, 58.0000000, 60.0000000, 62.0000000, 64.0000000, 66.0000000, 68.0000000, 70.0000000, 72.0000000, 74.0000000, 76.0000000, 78.0000000, 80.0000000, 82.0000000, 84.0000000, 86.0000000, 88.0000000, 90.0000000, 92.0000000, 94.0000000, 96.0000000, 98.0000000, 100.0000000, 102.0000000]

 (3)

NULL

N := numelems(x)

102

 (4)

``

`#mover(mi("U"),mo("&rarr;"))` := `<,>`(seq(U[n](t), n = 1 .. N))

_rtable[36893488152063361860]

 (5)

 

 

odes := seq(ode(n), n = 1 .. N)

ode(1), ode(2), ode(3), ode(4), ode(5), ode(6), ode(7), ode(8), ode(9), ode(10), ode(11), ode(12), ode(13), ode(14), ode(15), ode(16), ode(17), ode(18), ode(19), ode(20), ode(21), ode(22), ode(23), ode(24), ode(25), ode(26), ode(27), ode(28), ode(29), ode(30), ode(31), ode(32), ode(33), ode(34), ode(35), ode(36), ode(37), ode(38), ode(39), ode(40), ode(41), ode(42), ode(43), ode(44), ode(45), ode(46), ode(47), ode(48), ode(49), ode(50), ode(51), ode(52), ode(53), ode(54), ode(55), ode(56), ode(57), ode(58), ode(59), ode(60), ode(61), ode(62), ode(63), ode(64), ode(65), ode(66), ode(67), ode(68), ode(69), ode(70), ode(71), ode(72), ode(73), ode(74), ode(75), ode(76), ode(77), ode(78), ode(79), ode(80), ode(81), ode(82), ode(83), ode(84), ode(85), ode(86), ode(87), ode(88), ode(89), ode(90), ode(91), ode(92), ode(93), ode(94), ode(95), ode(96), ode(97), ode(98), ode(99), ode(100), ode(101), ode(102)

 (6)

 

W := proc (Un) options operator, arrow; -(1/4)*(b*d-e*q)^2/(b^2*(e+b*Un^2))-(2*b*c*m^2+(1/4)*q^2)*Un^2/b+c*Un^4 end proc; dW := D(W); dW(Un)

(1/2)*(b*d-e*q)^2*Un/(b*(Un^2*b+e)^2)-2*(2*b*c*m^2+(1/4)*q^2)*Un/b+4*c*Un^3

 (7)

#where each ode is defined by this equation:

"ode(n) := (U)[n] - `Delta__d`(n) + delta*(U)[n] = ((b d-e q)^2 U[n] )/(2 b (b U[n]^(2) +e)^2)-(2 (2 b c m^2+(q^2)/4) U[n] )/b+4 c U[n]^(3) "

proc (n) options operator, arrow, function_assign; (diff(`#mover(mi("U"),mo("&rarr;"))`, t, t))[n]-Delta__d(n)+`#mover(mi("&delta;",fontstyle = "normal"),mo("&circ;"))`*(diff(`#mover(mi("U"),mo("&rarr;"))`, t))[n] = (1/2)*(b*d-e*q)^2*`#mover(mi("U"),mo("&rarr;"))`[n]/(b*(b*`#mover(mi("U"),mo("&rarr;"))`[n]^2+e)^2)-(4*b*c*m^2+(1/2)*q^2)*`#mover(mi("U"),mo("&rarr;"))`[n]/b+4*c*`#mover(mi("U"),mo("&rarr;"))`[n]^3 end proc

 (8)

#and  the operator function is defined as:

"`Delta__d`(n) := {[[0,n=1 or n = N],[U[n+1] - 2*U[n] +U[n-1],otherwise]]; "

proc (n) options operator, arrow, function_assign; piecewise(n = 1 or n = N, 0, `#mover(mi("U"),mo("&rarr;"))`[n+1]-2*`#mover(mi("U"),mo("&rarr;"))`[n]+`#mover(mi("U"),mo("&rarr;"))`[n-1]) end proc

 (9)

 

#Define a sequence of initial conditions:

initPositions := seq(U[n](0) = U__0(n), n = 1 .. N); initVelocities := seq((D(U[n]))(0) = 0, n = 1 .. N)

(D(U[1]))(0) = 0, (D(U[2]))(0) = 0, (D(U[3]))(0) = 0, (D(U[4]))(0) = 0, (D(U[5]))(0) = 0, (D(U[6]))(0) = 0, (D(U[7]))(0) = 0, (D(U[8]))(0) = 0, (D(U[9]))(0) = 0, (D(U[10]))(0) = 0, (D(U[11]))(0) = 0, (D(U[12]))(0) = 0, (D(U[13]))(0) = 0, (D(U[14]))(0) = 0, (D(U[15]))(0) = 0, (D(U[16]))(0) = 0, (D(U[17]))(0) = 0, (D(U[18]))(0) = 0, (D(U[19]))(0) = 0, (D(U[20]))(0) = 0, (D(U[21]))(0) = 0, (D(U[22]))(0) = 0, (D(U[23]))(0) = 0, (D(U[24]))(0) = 0, (D(U[25]))(0) = 0, (D(U[26]))(0) = 0, (D(U[27]))(0) = 0, (D(U[28]))(0) = 0, (D(U[29]))(0) = 0, (D(U[30]))(0) = 0, (D(U[31]))(0) = 0, (D(U[32]))(0) = 0, (D(U[33]))(0) = 0, (D(U[34]))(0) = 0, (D(U[35]))(0) = 0, (D(U[36]))(0) = 0, (D(U[37]))(0) = 0, (D(U[38]))(0) = 0, (D(U[39]))(0) = 0, (D(U[40]))(0) = 0, (D(U[41]))(0) = 0, (D(U[42]))(0) = 0, (D(U[43]))(0) = 0, (D(U[44]))(0) = 0, (D(U[45]))(0) = 0, (D(U[46]))(0) = 0, (D(U[47]))(0) = 0, (D(U[48]))(0) = 0, (D(U[49]))(0) = 0, (D(U[50]))(0) = 0, (D(U[51]))(0) = 0, (D(U[52]))(0) = 0, (D(U[53]))(0) = 0, (D(U[54]))(0) = 0, (D(U[55]))(0) = 0, (D(U[56]))(0) = 0, (D(U[57]))(0) = 0, (D(U[58]))(0) = 0, (D(U[59]))(0) = 0, (D(U[60]))(0) = 0, (D(U[61]))(0) = 0, (D(U[62]))(0) = 0, (D(U[63]))(0) = 0, (D(U[64]))(0) = 0, (D(U[65]))(0) = 0, (D(U[66]))(0) = 0, (D(U[67]))(0) = 0, (D(U[68]))(0) = 0, (D(U[69]))(0) = 0, (D(U[70]))(0) = 0, (D(U[71]))(0) = 0, (D(U[72]))(0) = 0, (D(U[73]))(0) = 0, (D(U[74]))(0) = 0, (D(U[75]))(0) = 0, (D(U[76]))(0) = 0, (D(U[77]))(0) = 0, (D(U[78]))(0) = 0, (D(U[79]))(0) = 0, (D(U[80]))(0) = 0, (D(U[81]))(0) = 0, (D(U[82]))(0) = 0, (D(U[83]))(0) = 0, (D(U[84]))(0) = 0, (D(U[85]))(0) = 0, (D(U[86]))(0) = 0, (D(U[87]))(0) = 0, (D(U[88]))(0) = 0, (D(U[89]))(0) = 0, (D(U[90]))(0) = 0, (D(U[91]))(0) = 0, (D(U[92]))(0) = 0, (D(U[93]))(0) = 0, (D(U[94]))(0) = 0, (D(U[95]))(0) = 0, (D(U[96]))(0) = 0, (D(U[97]))(0) = 0, (D(U[98]))(0) = 0, (D(U[99]))(0) = 0, (D(U[100]))(0) = 0, (D(U[101]))(0) = 0, (D(U[102]))(0) = 0

 (10)

 

#where:

"`U__0`(n) :={[[0.0,n = 1 or n = N],[a*sech((h*x[n] ) ),otherwise]];"

proc (n) options operator, arrow, function_assign; piecewise(n = 1 or n = N, 0., a*sech(h*x[n])) end proc

 (11)

 

 

delta := .5; b := -0.222159366e-3; c := 1.088046084; d := 0.1308858509e-2; e := 0.394423507e-3; q := -0.64844084e-3; m := 1.518466566
 

a := 1.91; j := 1
 

NULLplot([seq([x[n], U__0(n)], n = 1 .. N)], style = point, axes = boxed, gridlines = true, labels = ['x__n', 'U__n(t = 0)'], labelfont = [Times, 16], size = [500, 300], color = "red")

  

sols := dsolve({odes, initPositions, initVelocities}, numeric, maxfun = 0, output = listprocedure)

`[Length of output exceeds limit of 1000000]`

 (12)

plots:-odeplot(sols, [t, U[67](t)], t = 0 .. .3000, color = "green")

  

`#mover(mi("Usol"),mo("&rarr;"))` := `<,>`(seq(eval(U[n](t), sols), n = 1 .. N))

`#mover(mi("Usol"),mo("&rarr;"))`[67](3.0)

Error, (in `#mover(mi("Usol"),mo("&rarr;"))`[67]) cannot evaluate the solution further right of .66258540, probably a singularity

  

" SnapshotPlot(t) := plot([seq([x[n], Usol[n](t)], n = 1..N)], style = point, symbolsize = 10, symbol = solidcircle, axes = boxed, gridlines = true, labels = ['`x__n`', '`U__n`'], labelfont = [Times, 15], title = typeset(`t = `, t), titlefont = [Garmond, 20], size = [600, 300],color="blue" ):"

Warning, (in SnapshotPlot) `n` is implicitly declared local

  

SnapshotPlot(0)

  

plots:-display(Array(1 .. 2, 1 .. 2, `<,>`(`<|>`(SnapshotPlot(0), SnapshotPlot(10)), `<|>`(SnapshotPlot(100), SnapshotPlot(3000)))))

Error, (in `#mover(mi("Usol"),mo("&rarr;"))`[1]) cannot evaluate the solution further right of .66258540, probably a singularity

  

NULL

NULL

 

 

 

 

NULL


 

Download v1.mw

@acer I probably deleted it when I saved the document last time. My bad

@acer I tried several values for ex...

Paras31 130
September 30 2024
0 0

@acer I tried several values for example `#mover(mi("&delta;",fontstyle = "normal"),mo("&circ;"))` := 0.1000000000.

and I had this result Error, (in `#mover(mi("Usol"),mo("&rarr;"))`[1]) cannot evaluate the solution further right of .62649992, probably a singularity

@acer Thank you very much for the p...

Paras31 130
September 26 2024
0 0

@acer Thank you very much for the presumptive answer. However, I was wondering if there is also in Maple a command that calculates the domain and the range of functions.

@janhardo If you run the following ...

Paras31 130
September 16 2024
0 0

@janhardo If you run the following worksheet for every amplitude I referred to above you will see that Maple gives the desired results. The Discrete Differential Equations are more difficult than ODEs. 

 

 I have created a procedure and I can plot the initial condition at t=0. Now I am trying to create a procedure to plot the results at t=3000 for different values of a.


 

The discrete Klein - Gordon equation (DKG) with friction is given by the equation:

 

diff(U(t), t, t)-K(`U__n+!`-2*U__n+`U__n-1`)+delta*U__n+(D(W))(U__n) = 0, beta > 0, delta > 0

In the above equation, W describes the pontential function:

W(U__n) = -(1/2)*`&omega;__d`^2*U__n^2+(1/4)*beta*`&omega;__d`^2*U__n^4

 

to which every coupled unit U__nadheres. In Eq. (1), the variable U__n(t)is the unknown displacement of the oscillator occupying the
n-th position of the lattice, and k = 1/h^2 is the discretization parameter. We denote by h the distance between the oscillators of the lattice. The chain (DKG) contains linear damping with a damping coefficient delta > 0, while β s the coefficient of the nonlinear cubic term.

U__n

 (1)


For the DKG chain (1), we will consider the problem of initial-boundary values, with initial conditions

 

U__n(0) = `U__n,0` and U__n(0) = `U__n,1` and `in`(`U__n,1`, 2*real^(k+2))

 

and Dirichlet boundary conditions at the boundary points x__0 = -(1/2)*Land "`x__k+1`=L/(2),"that is,

U__0 = `U__k+1` and `U__k+1` = 0, t >= 0*3

 

Therefore, when necessary, we will use the short notation `&Delta;__d`or the one-dimensional discrete Laplacian

`U__&Delta;__d__U__n&in;&Zopf;` = U__-2*U__n+4*`U__n-1`

 

We investigate numerically the dynamics of the system (1)-(2)-(3).Our first aim is to conduct a numerical study of the property of Dynamic Stability of the system,which directly depends on the existence and linear stability of the branches of equilibrium points. 

 

For the discussion of numerical results,it is also important to emphasize the role of the parameter `&omega;__d`By changing the time variable proc (t) options operator, arrow; 1/h end proc we rewrite Eq.(1) in the form

 

diff(U(t), t, t)-`&Delta;__d`*U(t)+(diff(delta(x), x))*U__n = `&Omega;__d`(-`&beta;U__n`^3+U__n)^2, t > 0, `&Omega;__d`^2 = h^2*`&omega;__d`^2, diff(delta(x), x) = `h&delta;`

 

The change in the scaling of the lattice parameter of the problem makes it important to comment here on the nature of the continuous and anti-continuous limit. For the anti-continuous limit, we need to consider the case in Eq. (1) where proc (`&omega;__d`) options operator, arrow; infinity end procIn the case of nonlinearity present in the governing equations, the continuous limit needs to be taken as well. On the other hand, for small values of the parameter, the system becomes significant. However, because of the model we consider, we take the asymptotic linear limit as

proc (`&omega;__d`) options operator, arrow; 0 end proc

proc (omega__d) options operator, arrow; 0 end proc

 (2)

 

 

We consider spatially extended initial conditions of the form:

U__n(0) = `U__n,0` and `U__n,0` = a*sin(`j&pi;hn`/L), j = 1, () .. K

 

where h = L/(K+1) is the distance of the grid and  a > 0 is the amplitude of the initial condition
We also assume zero initial velocity:

 

"(`U__n`(0)=`U__n,1`=0) "

 

 

> 

restart; K := 99; L := 200

 

Define spatial grid:

> 

x := [seq(-(1/2)*L+h*i, i = 0 .. K+1)]

Number of oscillators

> 

N := numelems(x)

 

`#mover(mi("U"),mo("&rarr;"))` is a vector of functions

> 

`#mover(mi("U"),mo("&rarr;"))` := `<,>`(seq(U[n](t), n = 1 .. N))

 

Create a sequence of odes

 

> 

odes := seq(ode(n), n = 1 .. N)

 

where each ode is defined by this equation:

> 

"ode(n) := (U)[n] - `Delta__d`(n) + delta*(U)[n] = `Omega__d`^(2)*(U[n] - beta*(U[n])^(3)): "

 

and  the operator function is defined as:

> 

"`Delta__d`(n) := {[[0,n=1 or n = N],[U[n+1] - 2*U[n] +U[n-1],otherwise]]:"

 

Define a sequence of initial conditions:

 

> 

initPositions := seq(U[n](0) = U__0(n), n = 1 .. N); initVelocities := seq((D(U[n]))(0) = 0, n = 1 .. N)

where:

> 

" `U__0`(n) :={[[0.0,n = 1 or n = N],[a*sin((j*Pi*h*x[n]/L) ),otherwise]]:"

 

However, I don't understand this initial condition function. Is j a value or one of several values?

 

Values of constants:

> 

beta := 1.0; delta := 0.5e-1; `&omega;__d` := sqrt(1); a := 1.91; j := 2

Some calculated constants:

> 

h := L/(K+1.0); `&Omega;__d` := h*`&omega;__d`; `#mover(mi("&delta;",fontstyle = "normal"),mo("&circ;"))` := h*delta

 

Plotting the initial position function: Is this correct?

 

> 

plot([seq([x[n], U__0(n)], n = 1 .. N)], style = point, axes = boxed, gridlines = true, labels = ['x__n', 'U__n(t = 0)'], labelfont = [Times, 16], size = [500, 300])

 

 

 

Solve but DO NOT DISPLAY.

 

> 

sols := dsolve({odes, initPositions, initVelocities}, numeric, maxfun = 0, output = listprocedure)

Plotting vs. time. Here the inefficiencies of the listprocedure option is noticable. Use odeplot

> 

plots:-odeplot(sols, [t, U[67](t)], t = 0 .. 200)

 

Extract out the solutions for convenience

> 

`#mover(mi("Usol"),mo("&rarr;"))` := `<,>`(seq(eval(U[n](t), sols), n = 1 .. N))

 

Looking a solution for oscillator 67 at t = 3.0

> 

`#mover(mi("Usol"),mo("&rarr;"))`[67](3.0)

HFloat(1.0740301674917803)

 (3)

Create the "snapshot in time" plot.

> 

"SnapshotPlot(t) := plot([seq([x[n], Usol[n](t)], n = 1..N)], style = point, symbolsize = 10, symbol = solidcircle, axes = boxed, gridlines = true, labels = ['`x__n`', '`U__n`'], labelfont = [Times, 15], title = typeset(`t = `, t), titlefont = [Garmond, 20], size = [600, 300]):"

Warning, (in SnapshotPlot) `n` is implicitly declared local

  

 

At t = 0, it matches the initial conditions.

> 

SnapshotPlot(0)

Sample Array plot

> 

plots:-display(Array(1 .. 2, 1 .. 2, `<,>`(`<|>`(SnapshotPlot(0), SnapshotPlot(10)), `<|>`(SnapshotPlot(50), SnapshotPlot(3000)))))

 

 
 

 

 

  

 

Some animation with time.

> 

plots:-animate(SnapshotPlot, [t], t = 0 .. 300, size = [600, 300], frames = 150, paraminfo = false)

 
> 

NULL
> 

NULL
> 

U__0 := proc (n, a) options operator, arrow; piecewise(n = 1 or n = N, 0., a*sin(2*j*Pi*x[n]/L)) end proc; initPositions := proc (a) local n; return seq(U__0(n, a), n = 1 .. N) end proc; initVelocities := seq((D(U[n]))(0) = 0, n = 1 .. N); a_values := [2.0, 1.95, 1.9, 1.85, 1.82, 1.81, 1.55, 1.46, 1.45]; colors := [red, blue, maroon, magenta, indigo, yellowgreen, black, gray, darkblue]; num_plots := nops(a_values); cols := 3; rows := ceil(num_plots/cols); plot_array := Array(1 .. rows, 1 .. cols); for i to num_plots do row := ceil(i/cols); col := i-(row-1)*cols; plot_array[row, col] := plot([seq([x[n], U__0(n, a_values[i])], n = 1 .. N)], style = point, axes = boxed, gridlines = true, labels = ['x__n', 'U__n(t = 0)'], labelfont = [Times, 16], title = cat("Plot for a = ", a_values[i]), size = [400, 200], color = colors[i]) end do; for r to rows do for c to cols do if `not`(assigned(plot_array[r, c])) then plot_array[r, c] := plot([]) end if end do end do; plots:-display(plot_array)

 

 
 

 

 

  
> 

NULL


 

Download DKG_EQUATION_V2.mw
 

The discrete Klein - Gordon equation (DKG) with friction is given by the equation:

 

diff(U(t), t, t)-K(`U__n+!`-2*U__n+`U__n-1`)+delta*U__n+(D(W))(U__n) = 0, beta > 0, delta > 0

In the above equation, W describes the pontential function:

W(U__n) = -(1/2)*`&omega;__d`^2*U__n^2+(1/4)*beta*`&omega;__d`^2*U__n^4

 

to which every coupled unit U__nadheres. In Eq. (1), the variable U__n(t)is the unknown displacement of the oscillator occupying the
n-th position of the lattice, and k = 1/h^2 is the discretization parameter. We denote by h the distance between the oscillators of the lattice. The chain (DKG) contains linear damping with a damping coefficient delta > 0, while β s the coefficient of the nonlinear cubic term.

U__n

 (1)


For the DKG chain (1), we will consider the problem of initial-boundary values, with initial conditions

 

U__n(0) = `U__n,0` and U__n(0) = `U__n,1` and `in`(`U__n,1`, 2*real^(k+2))

 

and Dirichlet boundary conditions at the boundary points x__0 = -(1/2)*Land "`x__k+1`=L/(2),"that is,

U__0 = `U__k+1` and `U__k+1` = 0, t >= 0*3

 

Therefore, when necessary, we will use the short notation `&Delta;__d`or the one-dimensional discrete Laplacian

`U__&Delta;__d__U__n&in;&Zopf;` = U__-2*U__n+4*`U__n-1`

 

We investigate numerically the dynamics of the system (1)-(2)-(3).Our first aim is to conduct a numerical study of the property of Dynamic Stability of the system,which directly depends on the existence and linear stability of the branches of equilibrium points. 

 

For the discussion of numerical results,it is also important to emphasize the role of the parameter `&omega;__d`By changing the time variable proc (t) options operator, arrow; 1/h end proc we rewrite Eq.(1) in the form

 

diff(U(t), t, t)-`&Delta;__d`*U(t)+(diff(delta(x), x))*U__n = `&Omega;__d`(-`&beta;U__n`^3+U__n)^2, t > 0, `&Omega;__d`^2 = h^2*`&omega;__d`^2, diff(delta(x), x) = `h&delta;`

 

The change in the scaling of the lattice parameter of the problem makes it important to comment here on the nature of the continuous and anti-continuous limit. For the anti-continuous limit, we need to consider the case in Eq. (1) where proc (`&omega;__d`) options operator, arrow; infinity end procIn the case of nonlinearity present in the governing equations, the continuous limit needs to be taken as well. On the other hand, for small values of the parameter, the system becomes significant. However, because of the model we consider, we take the asymptotic linear limit as

proc (`&omega;__d`) options operator, arrow; 0 end proc

proc (omega__d) options operator, arrow; 0 end proc

 (2)

 

 

We consider spatially extended initial conditions of the form:

U__n(0) = `U__n,0` and `U__n,0` = a*sin(`j&pi;hn`/L), j = 1, () .. K

 

where h = L/(K+1) is the distance of the grid and  a > 0 is the amplitude of the initial condition
We also assume zero initial velocity:

 

"(`U__n`(0)=`U__n,1`=0) "

 

 

> 

restart; K := 99; L := 200

 

Define spatial grid:

> 

x := [seq(-(1/2)*L+h*i, i = 0 .. K+1)]

Number of oscillators

> 

N := numelems(x)

 

`#mover(mi("U"),mo("&rarr;"))` is a vector of functions

> 

`#mover(mi("U"),mo("&rarr;"))` := `<,>`(seq(U[n](t), n = 1 .. N))

 

Create a sequence of odes

 

> 

odes := seq(ode(n), n = 1 .. N)

 

where each ode is defined by this equation:

> 

"ode(n) := (U)[n] - `Delta__d`(n) + delta*(U)[n] = `Omega__d`^(2)*(U[n] - beta*(U[n])^(3)): "

 

and  the operator function is defined as:

> 

"`Delta__d`(n) := {[[0,n=1 or n = N],[U[n+1] - 2*U[n] +U[n-1],otherwise]]:"

 

Define a sequence of initial conditions:

 

> 

initPositions := seq(U[n](0) = U__0(n), n = 1 .. N); initVelocities := seq((D(U[n]))(0) = 0, n = 1 .. N)

where:

> 

" `U__0`(n) :={[[0.0,n = 1 or n = N],[a*sin((j*Pi*h*x[n]/L) ),otherwise]]:"

 

However, I don't understand this initial condition function. Is j a value or one of several values?

 

Values of constants:

> 

beta := 1.0; delta := 0.5e-1; `&omega;__d` := sqrt(1); a := 1.91; j := 2

Some calculated constants:

> 

h := L/(K+1.0); `&Omega;__d` := h*`&omega;__d`; `#mover(mi("&delta;",fontstyle = "normal"),mo("&circ;"))` := h*delta

 

Plotting the initial position function: Is this correct?

 

> 

plot([seq([x[n], U__0(n)], n = 1 .. N)], style = point, axes = boxed, gridlines = true, labels = ['x__n', 'U__n(t = 0)'], labelfont = [Times, 16], size = [500, 300])

 

 

 

Solve but DO NOT DISPLAY.

 

> 

sols := dsolve({odes, initPositions, initVelocities}, numeric, maxfun = 0, output = listprocedure)

Plotting vs. time. Here the inefficiencies of the listprocedure option is noticable. Use odeplot

> 

plots:-odeplot(sols, [t, U[67](t)], t = 0 .. 200)

 

Extract out the solutions for convenience

> 

`#mover(mi("Usol"),mo("&rarr;"))` := `<,>`(seq(eval(U[n](t), sols), n = 1 .. N))

 

Looking a solution for oscillator 67 at t = 3.0

> 

`#mover(mi("Usol"),mo("&rarr;"))`[67](3.0)

HFloat(1.0740301674917803)

 (3)

Create the "snapshot in time" plot.

> 

"SnapshotPlot(t) := plot([seq([x[n], Usol[n](t)], n = 1..N)], style = point, symbolsize = 10, symbol = solidcircle, axes = boxed, gridlines = true, labels = ['`x__n`', '`U__n`'], labelfont = [Times, 15], title = typeset(`t = `, t), titlefont = [Garmond, 20], size = [600, 300]):"

Warning, (in SnapshotPlot) `n` is implicitly declared local

  

 

At t = 0, it matches the initial conditions.

> 

SnapshotPlot(0)

Sample Array plot

> 

plots:-display(Array(1 .. 2, 1 .. 2, `<,>`(`<|>`(SnapshotPlot(0), SnapshotPlot(10)), `<|>`(SnapshotPlot(50), SnapshotPlot(3000)))))

 

 
 

 

 

  

 

Some animation with time.

> 

plots:-animate(SnapshotPlot, [t], t = 0 .. 300, size = [600, 300], frames = 150, paraminfo = false)

 
> 

NULL
> 

NULL
> 

U__0 := proc (n, a) options operator, arrow; piecewise(n = 1 or n = N, 0., a*sin(2*j*Pi*x[n]/L)) end proc; initPositions := proc (a) local n; return seq(U__0(n, a), n = 1 .. N) end proc; initVelocities := seq((D(U[n]))(0) = 0, n = 1 .. N); a_values := [2.0, 1.95, 1.9, 1.85, 1.82, 1.81, 1.55, 1.46, 1.45]; colors := [red, blue, maroon, magenta, indigo, yellowgreen, black, gray, darkblue]; num_plots := nops(a_values); cols := 3; rows := ceil(num_plots/cols); plot_array := Array(1 .. rows, 1 .. cols); for i to num_plots do row := ceil(i/cols); col := i-(row-1)*cols; plot_array[row, col] := plot([seq([x[n], U__0(n, a_values[i])], n = 1 .. N)], style = point, axes = boxed, gridlines = true, labels = ['x__n', 'U__n(t = 0)'], labelfont = [Times, 16], title = cat("Plot for a = ", a_values[i]), size = [400, 200], color = colors[i]) end do; for r to rows do for c to cols do if `not`(assigned(plot_array[r, c])) then plot_array[r, c] := plot([]) end if end do end do; plots:-display(plot_array)

 

 
 

 

 

  
> 

NULL


 

Download DKG_EQUATION_V2.mw
 

 

 

 

@janhardo what do you mean?...

Paras31 130
September 15 2024
0 0

@janhardo what do you mean?

@janhardo I made some changes to ho...

Paras31 130
September 14 2024
0 0

@janhardo I made some changes to how to define the spatial grid and the initial position, Now I have the results I need. Now I have to figure out how to create a procedure to plot the solutions of different values of "a" at t=3000
 

Download DKG_EQUATION_V2.mw

@janhardo  After the help of P...

Paras31 130
September 14 2024
0 0

@janhardo  After the help of Professor @Scot Gould as shown in the worksheet we get the desired results. I agree it is not a simple system of ODEs

Paraskevopoulos problem

 

 

Counting constants:

> 

restart; K := 99; L := 200

 

Define spatial grid:

> 

x := evalf([seq(-(1/2)*L .. (1/2)*L, numelems = K+2)])

[-100., -98., -96., -94., -92., -90., -88., -86., -84., -82., -80., -78., -76., -74., -72., -70., -68., -66., -64., -62., -60., -58., -56., -54., -52., -50., -48., -46., -44., -42., -40., -38., -36., -34., -32., -30., -28., -26., -24., -22., -20., -18., -16., -14., -12., -10., -8., -6., -4., -2., 0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22., 24., 26., 28., 30., 32., 34., 36., 38., 40., 42., 44., 46., 48., 50., 52., 54., 56., 58., 60., 62., 64., 66., 68., 70., 72., 74., 76., 78., 80., 82., 84., 86., 88., 90., 92., 94., 96., 98., 100.]

 (1)

Number of oscillators

> 

N := numelems(x)

101

 (2)

 

`#mover(mi("U"),mo("&rarr;"))` is a vector of functions

> 

`#mover(mi("U"),mo("&rarr;"))` := `<,>`(seq(U[n](t), n = 1 .. N))

 

Create a sequence of odes

 

> 

odes := seq(ode(n), n = 1 .. N)

 

where each ode is defined by this equation:

> 

"ode(n) := (U)[n] - `Delta__d`(n) + delta*(U)[n] = `Omega__d`^(2)*(U[n] - beta*(U[n])^(3)): "

 

and  the operator function is defined as:

> 

"`Delta__d`(n) := {[[0,n=1 or n = N],[U[n+1] - 2*U[n] +U[n-1],otherwise]]:"

 

Define a sequence of initial conditions:

 

> 

initPositions := seq(U[n](0) = U__0(n), n = 1 .. N); initVelocities := seq((D(U[n]))(0) = 0, n = 1 .. N)

U[1](0) = 0., U[2](0) = 1.91*sin((2/25)*Pi), U[3](0) = 1.91*sin((3/25)*Pi), U[4](0) = 1.91*sin((4/25)*Pi), U[5](0) = 1.91*sin((1/5)*Pi), U[6](0) = 1.91*sin((6/25)*Pi), U[7](0) = 1.91*sin((7/25)*Pi), U[8](0) = 1.91*sin((8/25)*Pi), U[9](0) = 1.91*sin((9/25)*Pi), U[10](0) = 1.91*sin((2/5)*Pi), U[11](0) = 1.91*sin((11/25)*Pi), U[12](0) = 1.91*sin((12/25)*Pi), U[13](0) = 1.91*sin((12/25)*Pi), U[14](0) = 1.91*sin((11/25)*Pi), U[15](0) = 1.91*sin((2/5)*Pi), U[16](0) = 1.91*sin((9/25)*Pi), U[17](0) = 1.91*sin((8/25)*Pi), U[18](0) = 1.91*sin((7/25)*Pi), U[19](0) = 1.91*sin((6/25)*Pi), U[20](0) = 1.91*sin((1/5)*Pi), U[21](0) = 1.91*sin((4/25)*Pi), U[22](0) = 1.91*sin((3/25)*Pi), U[23](0) = 1.91*sin((2/25)*Pi), U[24](0) = 1.91*sin((1/25)*Pi), U[25](0) = 0., U[26](0) = -1.91*sin((1/25)*Pi), U[27](0) = -1.91*sin((2/25)*Pi), U[28](0) = -1.91*sin((3/25)*Pi), U[29](0) = -1.91*sin((4/25)*Pi), U[30](0) = -1.91*sin((1/5)*Pi), U[31](0) = -1.91*sin((6/25)*Pi), U[32](0) = -1.91*sin((7/25)*Pi), U[33](0) = -1.91*sin((8/25)*Pi), U[34](0) = -1.91*sin((9/25)*Pi), U[35](0) = -1.91*sin((2/5)*Pi), U[36](0) = -1.91*sin((11/25)*Pi), U[37](0) = -1.91*sin((12/25)*Pi), U[38](0) = -1.91*sin((12/25)*Pi), U[39](0) = -1.91*sin((11/25)*Pi), U[40](0) = -1.91*sin((2/5)*Pi), U[41](0) = -1.91*sin((9/25)*Pi), U[42](0) = -1.91*sin((8/25)*Pi), U[43](0) = -1.91*sin((7/25)*Pi), U[44](0) = -1.91*sin((6/25)*Pi), U[45](0) = -1.91*sin((1/5)*Pi), U[46](0) = -1.91*sin((4/25)*Pi), U[47](0) = -1.91*sin((3/25)*Pi), U[48](0) = -1.91*sin((2/25)*Pi), U[49](0) = -1.91*sin((1/25)*Pi), U[50](0) = 0., U[51](0) = 1.91*sin((1/25)*Pi), U[52](0) = 1.91*sin((2/25)*Pi), U[53](0) = 1.91*sin((3/25)*Pi), U[54](0) = 1.91*sin((4/25)*Pi), U[55](0) = 1.91*sin((1/5)*Pi), U[56](0) = 1.91*sin((6/25)*Pi), U[57](0) = 1.91*sin((7/25)*Pi), U[58](0) = 1.91*sin((8/25)*Pi), U[59](0) = 1.91*sin((9/25)*Pi), U[60](0) = 1.91*sin((2/5)*Pi), U[61](0) = 1.91*sin((11/25)*Pi), U[62](0) = 1.91*sin((12/25)*Pi), U[63](0) = 1.91*sin((12/25)*Pi), U[64](0) = 1.91*sin((11/25)*Pi), U[65](0) = 1.91*sin((2/5)*Pi), U[66](0) = 1.91*sin((9/25)*Pi), U[67](0) = 1.91*sin((8/25)*Pi), U[68](0) = 1.91*sin((7/25)*Pi), U[69](0) = 1.91*sin((6/25)*Pi), U[70](0) = 1.91*sin((1/5)*Pi), U[71](0) = 1.91*sin((4/25)*Pi), U[72](0) = 1.91*sin((3/25)*Pi), U[73](0) = 1.91*sin((2/25)*Pi), U[74](0) = 1.91*sin((1/25)*Pi), U[75](0) = 0., U[76](0) = -1.91*sin((1/25)*Pi), U[77](0) = -1.91*sin((2/25)*Pi), U[78](0) = -1.91*sin((3/25)*Pi), U[79](0) = -1.91*sin((4/25)*Pi), U[80](0) = -1.91*sin((1/5)*Pi), U[81](0) = -1.91*sin((6/25)*Pi), U[82](0) = -1.91*sin((7/25)*Pi), U[83](0) = -1.91*sin((8/25)*Pi), U[84](0) = -1.91*sin((9/25)*Pi), U[85](0) = -1.91*sin((2/5)*Pi), U[86](0) = -1.91*sin((11/25)*Pi), U[87](0) = -1.91*sin((12/25)*Pi), U[88](0) = -1.91*sin((12/25)*Pi), U[89](0) = -1.91*sin((11/25)*Pi), U[90](0) = -1.91*sin((2/5)*Pi), U[91](0) = -1.91*sin((9/25)*Pi), U[92](0) = -1.91*sin((8/25)*Pi), U[93](0) = -1.91*sin((7/25)*Pi), U[94](0) = -1.91*sin((6/25)*Pi), U[95](0) = -1.91*sin((1/5)*Pi), U[96](0) = -1.91*sin((4/25)*Pi), U[97](0) = -1.91*sin((3/25)*Pi), U[98](0) = -1.91*sin((2/25)*Pi), U[99](0) = -1.91*sin((1/25)*Pi), U[100](0) = 0., U[101](0) = 0.

 (3)

where:

> 

" `U__0`(n) :={[[0.0,n = 1 or n = N],[a*sin((j*Pi*2* n)/L),otherwise]]"

proc (n) options operator, arrow, function_assign; piecewise(n = 1 or n = N, 0., a*sin(2*j*Pi*n/L)) end proc

 (4)

 

However, I don't understand this initial condition function. Is j a value or one of several values?

 

Values of constants:

> 

beta := 1.0; delta := 0.5e-1; `&omega;__d` := sqrt(1); a := 1.91; j := 4

Some calculated constants:

> 

h := L/(K+1.0); `&Omega;__d` := h*`&omega;__d`; `#mover(mi("&delta;",fontstyle = "normal"),mo("&circ;"))` := h*delta

 

Plotting the initial position function: Is this correct?

 

> 

plot([seq([x[n], U__0(n)], n = 1 .. N)], style = point, axes = boxed, gridlines = true, labels = ['x__n', 'U__n(t = 0)'], labelfont = [Times, 16], size = [500, 300])

 

U__0(2)

1.91*sin((2/25)*Pi)

 (5)

 

Solve but DO NOT DISPLAY.

 

> 

sols := dsolve({odes, initPositions, initVelocities}, numeric, maxfun = 0, output = listprocedure)

Plotting vs. time. Here the inefficiencies of the listprocedure option is noticable. Use odeplot

> 

plots:-odeplot(sols, [t, U[67](t)], t = 0 .. 100)

 

Extract out the solutions for convenience

> 

`#mover(mi("Usol"),mo("&rarr;"))` := `<,>`(seq(eval(U[n](t), sols), n = 1 .. N))

 

Looking a solution for oscillator 67 at t = 3.0

> 

`#mover(mi("Usol"),mo("&rarr;"))`[67](3.0)

HFloat(0.20725065695143008)

 (6)

Create the "snapshot in time" plot.

> 

"SnapshotPlot(t) := plot([seq([x[n], Usol[n](t)], n = 1..N)], style = point, symbolsize = 10, symbol = solidcircle, axes = boxed, gridlines = true, labels = ['`x__n`', '`U__n`'], labelfont = [Times, 15], title = typeset(`t = `, t), titlefont = [Garmond, 20], size = [600, 300]):"

Warning, (in SnapshotPlot) `n` is implicitly declared local

  

 

At t = 0, it matches the initial conditions.

> 

SnapshotPlot(0)

 

Sample Array plot

> 

plots:-display(Array(1 .. 2, 1 .. 2, `<,>`(`<|>`(SnapshotPlot(0), SnapshotPlot(10)), `<|>`(SnapshotPlot(50), SnapshotPlot(250)))))

 

 
 

 

 

  

 

Some animation with time.

> 

plots:-animate(SnapshotPlot, [t], t = 0 .. 100, size = [600, 300], frames = 50, paraminfo = false)

 
> 

``
> 

NULL
> 

NULL
> 

NULL


 

Download DKG_Equation.mw
 

 

@janhardo To be more specific, it i...

Paras31 130
September 12 2024
0 0

@janhardo To be more specific, it is the Discrete Klein-Gordon Equation. It is a second-order differential equation. I have also worked on Matlab and Mathematica. But on Maple I face difficulties

@Hullzie16  I tried now to replicat...

Paras31 130
September 11 2024
0 0

@Hullzie16  I tried now to replicate it with a loop but this way failed too,

> 

restart;
with(plots):
> 

L := 200;
K := 99;
kappa := 1;
omegaD := 1;
beta := 1;
delta := 0.05;
j := 2;
dt:=0.05:
tmax := 3000;
h := L/(K + 1);
nsp := [(-L/2 + h*i) $ (i = 0 .. L/h)]:
km := nops(nsp);
omegaD2 := h^2*omegaD^2;
deltaHat := h*delta;
a := 1.8;
var := seq(x[i](t), i = 1 .. km):
initialPositions := seq(x[i](0) = a*sin(j*h*Pi*nsp[i]/L), i = 1 .. km);
initialVelocities := seq(D(x[i])(0) = 0, i = 1 .. km):

200

  

99

  

1

  

1

  

1

  

0.5e-1

  

2

  

3000

  

2

  

101

  

4

  

.10

  

1.8

  

x[1](0) = 0., x[2](0) = 1.8*sin((1/25)*Pi), x[3](0) = 1.8*sin((2/25)*Pi), x[4](0) = 1.8*sin((3/25)*Pi), x[5](0) = 1.8*sin((4/25)*Pi), x[6](0) = 1.8*sin((1/5)*Pi), x[7](0) = 1.8*sin((6/25)*Pi), x[8](0) = 1.8*sin((7/25)*Pi), x[9](0) = 1.8*sin((8/25)*Pi), x[10](0) = 1.8*sin((9/25)*Pi), x[11](0) = 1.8*sin((2/5)*Pi), x[12](0) = 1.8*sin((11/25)*Pi), x[13](0) = 1.8*sin((12/25)*Pi), x[14](0) = 1.8*sin((12/25)*Pi), x[15](0) = 1.8*sin((11/25)*Pi), x[16](0) = 1.8*sin((2/5)*Pi), x[17](0) = 1.8*sin((9/25)*Pi), x[18](0) = 1.8*sin((8/25)*Pi), x[19](0) = 1.8*sin((7/25)*Pi), x[20](0) = 1.8*sin((6/25)*Pi), x[21](0) = 1.8*sin((1/5)*Pi), x[22](0) = 1.8*sin((4/25)*Pi), x[23](0) = 1.8*sin((3/25)*Pi), x[24](0) = 1.8*sin((2/25)*Pi), x[25](0) = 1.8*sin((1/25)*Pi), x[26](0) = 0., x[27](0) = -1.8*sin((1/25)*Pi), x[28](0) = -1.8*sin((2/25)*Pi), x[29](0) = -1.8*sin((3/25)*Pi), x[30](0) = -1.8*sin((4/25)*Pi), x[31](0) = -1.8*sin((1/5)*Pi), x[32](0) = -1.8*sin((6/25)*Pi), x[33](0) = -1.8*sin((7/25)*Pi), x[34](0) = -1.8*sin((8/25)*Pi), x[35](0) = -1.8*sin((9/25)*Pi), x[36](0) = -1.8*sin((2/5)*Pi), x[37](0) = -1.8*sin((11/25)*Pi), x[38](0) = -1.8*sin((12/25)*Pi), x[39](0) = -1.8*sin((12/25)*Pi), x[40](0) = -1.8*sin((11/25)*Pi), x[41](0) = -1.8*sin((2/5)*Pi), x[42](0) = -1.8*sin((9/25)*Pi), x[43](0) = -1.8*sin((8/25)*Pi), x[44](0) = -1.8*sin((7/25)*Pi), x[45](0) = -1.8*sin((6/25)*Pi), x[46](0) = -1.8*sin((1/5)*Pi), x[47](0) = -1.8*sin((4/25)*Pi), x[48](0) = -1.8*sin((3/25)*Pi), x[49](0) = -1.8*sin((2/25)*Pi), x[50](0) = -1.8*sin((1/25)*Pi), x[51](0) = 0., x[52](0) = 1.8*sin((1/25)*Pi), x[53](0) = 1.8*sin((2/25)*Pi), x[54](0) = 1.8*sin((3/25)*Pi), x[55](0) = 1.8*sin((4/25)*Pi), x[56](0) = 1.8*sin((1/5)*Pi), x[57](0) = 1.8*sin((6/25)*Pi), x[58](0) = 1.8*sin((7/25)*Pi), x[59](0) = 1.8*sin((8/25)*Pi), x[60](0) = 1.8*sin((9/25)*Pi), x[61](0) = 1.8*sin((2/5)*Pi), x[62](0) = 1.8*sin((11/25)*Pi), x[63](0) = 1.8*sin((12/25)*Pi), x[64](0) = 1.8*sin((12/25)*Pi), x[65](0) = 1.8*sin((11/25)*Pi), x[66](0) = 1.8*sin((2/5)*Pi), x[67](0) = 1.8*sin((9/25)*Pi), x[68](0) = 1.8*sin((8/25)*Pi), x[69](0) = 1.8*sin((7/25)*Pi), x[70](0) = 1.8*sin((6/25)*Pi), x[71](0) = 1.8*sin((1/5)*Pi), x[72](0) = 1.8*sin((4/25)*Pi), x[73](0) = 1.8*sin((3/25)*Pi), x[74](0) = 1.8*sin((2/25)*Pi), x[75](0) = 1.8*sin((1/25)*Pi), x[76](0) = 0., x[77](0) = -1.8*sin((1/25)*Pi), x[78](0) = -1.8*sin((2/25)*Pi), x[79](0) = -1.8*sin((3/25)*Pi), x[80](0) = -1.8*sin((4/25)*Pi), x[81](0) = -1.8*sin((1/5)*Pi), x[82](0) = -1.8*sin((6/25)*Pi), x[83](0) = -1.8*sin((7/25)*Pi), x[84](0) = -1.8*sin((8/25)*Pi), x[85](0) = -1.8*sin((9/25)*Pi), x[86](0) = -1.8*sin((2/5)*Pi), x[87](0) = -1.8*sin((11/25)*Pi), x[88](0) = -1.8*sin((12/25)*Pi), x[89](0) = -1.8*sin((12/25)*Pi), x[90](0) = -1.8*sin((11/25)*Pi), x[91](0) = -1.8*sin((2/5)*Pi), x[92](0) = -1.8*sin((9/25)*Pi), x[93](0) = -1.8*sin((8/25)*Pi), x[94](0) = -1.8*sin((7/25)*Pi), x[95](0) = -1.8*sin((6/25)*Pi), x[96](0) = -1.8*sin((1/5)*Pi), x[97](0) = -1.8*sin((4/25)*Pi), x[98](0) = -1.8*sin((3/25)*Pi), x[99](0) = -1.8*sin((2/25)*Pi), x[100](0) = -1.8*sin((1/25)*Pi), x[101](0) = 0.

 (1)
> 

boundaryConditions := [x[1](t) = 0, x[km](t) = 0]

[x[1](t) = 0, x[101](t) = 0]

 (2)
> 

equations := seq(diff(x[n](t), t, t) - kappa*(x[n + 1](t) - 2*x[n](t) + x[n - 1](t)) + deltaHat*diff(x[n](t), t) - omegaD2^2*(x[n](t) - beta*x[n](t)^3) = 0, n = 2 .. km - 1), diff(x[1](t), t, t) = 0, diff(x[101](t), t, t) = 0:
> 

sol := dsolve({equations, initialPositions, initialVelocities}, maxfun=0,numeric, method = rkf45,output=listprocedure)

`[Length of output exceeds limit of 1000000]`

 (3)
> 

# Use odeplot to plot the trajectory of x[67](t)
odeplot(sol, [[t, x[67](t)]], t = 0 .. 100, numpoints = 3500, title = "Trajectory of x[67](t)", axes = boxed, gridlines=true, color = ["Blue"], size = [500, 300]);

 
> 

# Plot x[67](t) vs u[67](t) = diff(x[67](t), t)
odeplot(sol, [[x[67](t), diff(x[67](t), t)]], t = 0 .. 100, numpoints = 3500, title = "x_n vs u_n", axes = boxed, gridlines=true, color = ["Blue"], size = [500, 300]);

 
> 

odeplot(sol, [[t, diff(x[67](t), t)]], t = 0 .. 100, numpoints = 3500, title = "u_n", axes = boxed, gridlines = true, color = ["Red"], size = [500, 300])

 
> 

getDisplacements := proc (t_val) local U, i; U := [seq(eval(x[i](t_val), sol), i = 1 .. km)]; return U end proc
> 

TT := table(); for t_val from 0 by dt to 3000 do U := getDisplacements(t_val); plotU := pointplot([seq([nsp[i], U[i]], i = 1 .. K+2)], style = point, symbol = solidcircle, color = blue, title = cat("t = ", t_val, ", a = ", a), labels = ["n", "Displacement"], grid = [true, true], view = [-100 .. 100, -2 .. 2]); TT[t_val] := plotU; DocumentTools:-Tabulate([plotU]) end do

table( [ ] )

  

[x[1](0), x[2](0), x[3](0), x[4](0), x[5](0), x[6](0), x[7](0), x[8](0), x[9](0), x[10](0), x[11](0), x[12](0), x[13](0), x[14](0), x[15](0), x[16](0), x[17](0), x[18](0), x[19](0), x[20](0), x[21](0), x[22](0), x[23](0), x[24](0), x[25](0), x[26](0), x[27](0), x[28](0), x[29](0), x[30](0), x[31](0), x[32](0), x[33](0), x[34](0), x[35](0), x[36](0), x[37](0), x[38](0), x[39](0), x[40](0), x[41](0), x[42](0), x[43](0), x[44](0), x[45](0), x[46](0), x[47](0), x[48](0), x[49](0), x[50](0), x[51](0), x[52](0), x[53](0), x[54](0), x[55](0), x[56](0), x[57](0), x[58](0), x[59](0), x[60](0), x[61](0), x[62](0), x[63](0), x[64](0), x[65](0), x[66](0), x[67](0), x[68](0), x[69](0), x[70](0), x[71](0), x[72](0), x[73](0), x[74](0), x[75](0), x[76](0), x[77](0), x[78](0), x[79](0), x[80](0), x[81](0), x[82](0), x[83](0), x[84](0), x[85](0), x[86](0), x[87](0), x[88](0), x[89](0), x[90](0), x[91](0), x[92](0), x[93](0), x[94](0), x[95](0), x[96](0), x[97](0), x[98](0), x[99](0), x[100](0), x[101](0)]

  

Error, (in plots:-pointplot) points cannot be converted to floating-point values

  
 

NULL

Download v4.mw

pointplot need it...

Paras31 130
September 11 2024
0 0

@Hullzie16 As I can understand I have to create a pointplot spatial points over var:= seq(x[i](t), i = 1 .. km) the, but still I have not figured out how to replicate it

> 

restart;
with(plots):
> 

L := 200;
K := 99;
kappa := 1;
omegaD := 1;
beta := 1;
delta := 0.05;
j := 2;
tmax := 3000;
h := L/(K + 1);
nsp := [(-L/2 + h*i) $ (i = 0 .. L/h)]:
km := nops(nsp);
omegaD2 := h^2*omegaD^2;
deltaHat := h*delta;
a := 1.8;
var := seq(x[i](t), i = 1 .. km):
initialPositions := seq(x[i](0) = a*sin(j*h*Pi*nsp[i]/L), i = 1 .. km):
initialVelocities := seq(D(x[i])(0) = 0, i = 1 .. km):

200

  

99

  

1

  

1

  

1

  

0.5e-1

  

2

  

3000

  

2

  

101

  

4

  

.10

  

1.8

 (1)
> 

boundaryConditions := [x[1](t) = 0, x[km](t) = 0]

[x[1](t) = 0, x[101](t) = 0]

 (2)
> 

equations := seq(diff(x[n](t), t, t) - kappa*(x[n + 1](t) - 2*x[n](t) + x[n - 1](t)) + deltaHat*diff(x[n](t), t) - omegaD2^2*(x[n](t) - beta*x[n](t)^3) = 0, n = 2 .. km - 1), diff(x[1](t), t, t) = 0, diff(x[101](t), t, t) = 0:
> 

sol := dsolve({equations, initialPositions, initialVelocities}, maxfun=0,numeric, method = rkf45,output=listprocedure)

`[Length of output exceeds limit of 1000000]`

 (3)
> 

# Use odeplot to plot the trajectory of x[67](t)
odeplot(sol, [[t, x[67](t)]], t = 0 .. 100, numpoints = 3500, title = "Trajectory of x[67](t)", axes = boxed, gridlines=true, color = ["Blue"], size = [500, 300]);

 
> 

# Plot x[67](t) vs u[67](t) = diff(x[67](t), t)
odeplot(sol, [[x[67](t), diff(x[67](t), t)]], t = 0 .. 100, numpoints = 3500, title = "x_n vs u_n", axes = boxed, gridlines=true, color = ["Blue"], size = [500, 300]);

 
> 

odeplot(sol, [[t, diff(x[67](t), t)]], t = 0 .. 100, numpoints = 3500, title = "u_n", axes = boxed, gridlines = true, color = ["Red"], size = [500, 300])

 
> 

solVals := sol(0); solVals; dataPoints := [seq([nsp[i], evalf(solVals[i][2])], i = 1 .. km)]

[t(0) = HFloat(0.0), (x[1](t))(0) = HFloat(0.0), (diff(x[1](t), t))(0) = HFloat(0.0), (x[2](t))(0) = HFloat(0.225599820415749), (diff(x[2](t), t))(0) = HFloat(0.0), (x[3](t))(0) = HFloat(0.447641796896737), (diff(x[3](t), t))(0) = HFloat(0.0), (x[4](t))(0) = HFloat(0.66262419483242), (diff(x[4](t), t))(0) = HFloat(0.0), (x[5](t))(0) = HFloat(0.867156613383085), (diff(x[5](t), t))(0) = HFloat(0.0), (x[6](t))(0) = HFloat(1.05801345412645), (diff(x[6](t), t))(0) = HFloat(0.0), (x[7](t))(0) = HFloat(1.23218479067164), (diff(x[7](t), t))(0) = HFloat(0.0), (x[8](t))(0) = HFloat(1.38692383699642), (diff(x[8](t), t))(0) = HFloat(0.0), (x[9](t))(0) = HFloat(1.51979026590362), (diff(x[9](t), t))(0) = HFloat(0.0), (x[10](t))(0) = HFloat(1.62868869443883), (diff(x[10](t), t))(0) = HFloat(0.0), (x[11](t))(0) = HFloat(1.71190172933128), (diff(x[11](t), t))(0) = HFloat(0.0), (x[12](t))(0) = HFloat(1.76811705131164), (diff(x[12](t), t))(0) = HFloat(0.0), (x[13](t))(0) = HFloat(1.79644811117089), (diff(x[13](t), t))(0) = HFloat(0.0), (x[14](t))(0) = HFloat(1.79644811117089), (diff(x[14](t), t))(0) = HFloat(0.0), (x[15](t))(0) = HFloat(1.76811705131164), (diff(x[15](t), t))(0) = HFloat(0.0), (x[16](t))(0) = HFloat(1.71190172933128), (diff(x[16](t), t))(0) = HFloat(0.0), (x[17](t))(0) = HFloat(1.62868869443883), (diff(x[17](t), t))(0) = HFloat(0.0), (x[18](t))(0) = HFloat(1.51979026590362), (diff(x[18](t), t))(0) = HFloat(0.0), (x[19](t))(0) = HFloat(1.38692383699642), (diff(x[19](t), t))(0) = HFloat(0.0), (x[20](t))(0) = HFloat(1.23218479067164), (diff(x[20](t), t))(0) = HFloat(0.0), (x[21](t))(0) = HFloat(1.05801345412645), (diff(x[21](t), t))(0) = HFloat(0.0), (x[22](t))(0) = HFloat(0.867156613383085), (diff(x[22](t), t))(0) = HFloat(0.0), (x[23](t))(0) = HFloat(0.66262419483242), (diff(x[23](t), t))(0) = HFloat(0.0), (x[24](t))(0) = HFloat(0.447641796896737), (diff(x[24](t), t))(0) = HFloat(0.0), (x[25](t))(0) = HFloat(0.225599820415749), (diff(x[25](t), t))(0) = HFloat(0.0), (x[26](t))(0) = HFloat(0.0), (diff(x[26](t), t))(0) = HFloat(0.0), (x[27](t))(0) = HFloat(-0.225599820415749), (diff(x[27](t), t))(0) = HFloat(0.0), (x[28](t))(0) = HFloat(-0.447641796896737), (diff(x[28](t), t))(0) = HFloat(0.0), (x[29](t))(0) = HFloat(-0.66262419483242), (diff(x[29](t), t))(0) = HFloat(0.0), (x[30](t))(0) = HFloat(-0.867156613383085), (diff(x[30](t), t))(0) = HFloat(0.0), (x[31](t))(0) = HFloat(-1.05801345412645), (diff(x[31](t), t))(0) = HFloat(0.0), (x[32](t))(0) = HFloat(-1.23218479067164), (diff(x[32](t), t))(0) = HFloat(0.0), (x[33](t))(0) = HFloat(-1.38692383699642), (diff(x[33](t), t))(0) = HFloat(0.0), (x[34](t))(0) = HFloat(-1.51979026590362), (diff(x[34](t), t))(0) = HFloat(0.0), (x[35](t))(0) = HFloat(-1.62868869443883), (diff(x[35](t), t))(0) = HFloat(0.0), (x[36](t))(0) = HFloat(-1.71190172933128), (diff(x[36](t), t))(0) = HFloat(0.0), (x[37](t))(0) = HFloat(-1.76811705131164), (diff(x[37](t), t))(0) = HFloat(0.0), (x[38](t))(0) = HFloat(-1.79644811117089), (diff(x[38](t), t))(0) = HFloat(0.0), (x[39](t))(0) = HFloat(-1.79644811117089), (diff(x[39](t), t))(0) = HFloat(0.0), (x[40](t))(0) = HFloat(-1.76811705131164), (diff(x[40](t), t))(0) = HFloat(0.0), (x[41](t))(0) = HFloat(-1.71190172933128), (diff(x[41](t), t))(0) = HFloat(0.0), (x[42](t))(0) = HFloat(-1.62868869443883), (diff(x[42](t), t))(0) = HFloat(0.0), (x[43](t))(0) = HFloat(-1.51979026590362), (diff(x[43](t), t))(0) = HFloat(0.0), (x[44](t))(0) = HFloat(-1.38692383699642), (diff(x[44](t), t))(0) = HFloat(0.0), (x[45](t))(0) = HFloat(-1.23218479067164), (diff(x[45](t), t))(0) = HFloat(0.0), (x[46](t))(0) = HFloat(-1.05801345412645), (diff(x[46](t), t))(0) = HFloat(0.0), (x[47](t))(0) = HFloat(-0.867156613383085), (diff(x[47](t), t))(0) = HFloat(0.0), (x[48](t))(0) = HFloat(-0.66262419483242), (diff(x[48](t), t))(0) = HFloat(0.0), (x[49](t))(0) = HFloat(-0.447641796896737), (diff(x[49](t), t))(0) = HFloat(0.0), (x[50](t))(0) = HFloat(-0.225599820415749), (diff(x[50](t), t))(0) = HFloat(0.0), (x[51](t))(0) = HFloat(0.0), (diff(x[51](t), t))(0) = HFloat(0.0), (x[52](t))(0) = HFloat(0.225599820415749), (diff(x[52](t), t))(0) = HFloat(0.0), (x[53](t))(0) = HFloat(0.447641796896737), (diff(x[53](t), t))(0) = HFloat(0.0), (x[54](t))(0) = HFloat(0.66262419483242), (diff(x[54](t), t))(0) = HFloat(0.0), (x[55](t))(0) = HFloat(0.867156613383085), (diff(x[55](t), t))(0) = HFloat(0.0), (x[56](t))(0) = HFloat(1.05801345412645), (diff(x[56](t), t))(0) = HFloat(0.0), (x[57](t))(0) = HFloat(1.23218479067164), (diff(x[57](t), t))(0) = HFloat(0.0), (x[58](t))(0) = HFloat(1.38692383699642), (diff(x[58](t), t))(0) = HFloat(0.0), (x[59](t))(0) = HFloat(1.51979026590362), (diff(x[59](t), t))(0) = HFloat(0.0), (x[60](t))(0) = HFloat(1.62868869443883), (diff(x[60](t), t))(0) = HFloat(0.0), (x[61](t))(0) = HFloat(1.71190172933128), (diff(x[61](t), t))(0) = HFloat(0.0), (x[62](t))(0) = HFloat(1.76811705131164), (diff(x[62](t), t))(0) = HFloat(0.0), (x[63](t))(0) = HFloat(1.79644811117089), (diff(x[63](t), t))(0) = HFloat(0.0), (x[64](t))(0) = HFloat(1.79644811117089), (diff(x[64](t), t))(0) = HFloat(0.0), (x[65](t))(0) = HFloat(1.76811705131164), (diff(x[65](t), t))(0) = HFloat(0.0), (x[66](t))(0) = HFloat(1.71190172933128), (diff(x[66](t), t))(0) = HFloat(0.0), (x[67](t))(0) = HFloat(1.62868869443883), (diff(x[67](t), t))(0) = HFloat(0.0), (x[68](t))(0) = HFloat(1.51979026590362), (diff(x[68](t), t))(0) = HFloat(0.0), (x[69](t))(0) = HFloat(1.38692383699642), (diff(x[69](t), t))(0) = HFloat(0.0), (x[70](t))(0) = HFloat(1.23218479067164), (diff(x[70](t), t))(0) = HFloat(0.0), (x[71](t))(0) = HFloat(1.05801345412645), (diff(x[71](t), t))(0) = HFloat(0.0), (x[72](t))(0) = HFloat(0.867156613383085), (diff(x[72](t), t))(0) = HFloat(0.0), (x[73](t))(0) = HFloat(0.66262419483242), (diff(x[73](t), t))(0) = HFloat(0.0), (x[74](t))(0) = HFloat(0.447641796896737), (diff(x[74](t), t))(0) = HFloat(0.0), (x[75](t))(0) = HFloat(0.225599820415749), (diff(x[75](t), t))(0) = HFloat(0.0), (x[76](t))(0) = HFloat(0.0), (diff(x[76](t), t))(0) = HFloat(0.0), (x[77](t))(0) = HFloat(-0.225599820415749), (diff(x[77](t), t))(0) = HFloat(0.0), (x[78](t))(0) = HFloat(-0.447641796896737), (diff(x[78](t), t))(0) = HFloat(0.0), (x[79](t))(0) = HFloat(-0.66262419483242), (diff(x[79](t), t))(0) = HFloat(0.0), (x[80](t))(0) = HFloat(-0.867156613383085), (diff(x[80](t), t))(0) = HFloat(0.0), (x[81](t))(0) = HFloat(-1.05801345412645), (diff(x[81](t), t))(0) = HFloat(0.0), (x[82](t))(0) = HFloat(-1.23218479067164), (diff(x[82](t), t))(0) = HFloat(0.0), (x[83](t))(0) = HFloat(-1.38692383699642), (diff(x[83](t), t))(0) = HFloat(0.0), (x[84](t))(0) = HFloat(-1.51979026590362), (diff(x[84](t), t))(0) = HFloat(0.0), (x[85](t))(0) = HFloat(-1.62868869443883), (diff(x[85](t), t))(0) = HFloat(0.0), (x[86](t))(0) = HFloat(-1.71190172933128), (diff(x[86](t), t))(0) = HFloat(0.0), (x[87](t))(0) = HFloat(-1.76811705131164), (diff(x[87](t), t))(0) = HFloat(0.0), (x[88](t))(0) = HFloat(-1.79644811117089), (diff(x[88](t), t))(0) = HFloat(0.0), (x[89](t))(0) = HFloat(-1.79644811117089), (diff(x[89](t), t))(0) = HFloat(0.0), (x[90](t))(0) = HFloat(-1.76811705131164), (diff(x[90](t), t))(0) = HFloat(0.0), (x[91](t))(0) = HFloat(-1.71190172933128), (diff(x[91](t), t))(0) = HFloat(0.0), (x[92](t))(0) = HFloat(-1.62868869443883), (diff(x[92](t), t))(0) = HFloat(0.0), (x[93](t))(0) = HFloat(-1.51979026590362), (diff(x[93](t), t))(0) = HFloat(0.0), (x[94](t))(0) = HFloat(-1.38692383699642), (diff(x[94](t), t))(0) = HFloat(0.0), (x[95](t))(0) = HFloat(-1.23218479067164), (diff(x[95](t), t))(0) = HFloat(0.0), (x[96](t))(0) = HFloat(-1.05801345412645), (diff(x[96](t), t))(0) = HFloat(0.0), (x[97](t))(0) = HFloat(-0.867156613383085), (diff(x[97](t), t))(0) = HFloat(0.0), (x[98](t))(0) = HFloat(-0.66262419483242), (diff(x[98](t), t))(0) = HFloat(0.0), (x[99](t))(0) = HFloat(-0.447641796896737), (diff(x[99](t), t))(0) = HFloat(0.0), (x[100](t))(0) = HFloat(-0.225599820415749), (diff(x[100](t), t))(0) = HFloat(0.0), (x[101](t))(0) = HFloat(0.0), (diff(x[101](t), t))(0) = HFloat(0.0)]

  

[[-100, (t(0) = HFloat(0.0))[2]], [-98, ((x[1](t))(0) = HFloat(0.0))[2]], [-96, ((diff(x[1](t), t))(0) = HFloat(0.0))[2]], [-94, ((x[2](t))(0) = HFloat(0.225599820415749))[2]], [-92, ((diff(x[2](t), t))(0) = HFloat(0.0))[2]], [-90, ((x[3](t))(0) = HFloat(0.447641796896737))[2]], [-88, ((diff(x[3](t), t))(0) = HFloat(0.0))[2]], [-86, ((x[4](t))(0) = HFloat(0.66262419483242))[2]], [-84, ((diff(x[4](t), t))(0) = HFloat(0.0))[2]], [-82, ((x[5](t))(0) = HFloat(0.867156613383085))[2]], [-80, ((diff(x[5](t), t))(0) = HFloat(0.0))[2]], [-78, ((x[6](t))(0) = HFloat(1.05801345412645))[2]], [-76, ((diff(x[6](t), t))(0) = HFloat(0.0))[2]], [-74, ((x[7](t))(0) = HFloat(1.23218479067164))[2]], [-72, ((diff(x[7](t), t))(0) = HFloat(0.0))[2]], [-70, ((x[8](t))(0) = HFloat(1.38692383699642))[2]], [-68, ((diff(x[8](t), t))(0) = HFloat(0.0))[2]], [-66, ((x[9](t))(0) = HFloat(1.51979026590362))[2]], [-64, ((diff(x[9](t), t))(0) = HFloat(0.0))[2]], [-62, ((x[10](t))(0) = HFloat(1.62868869443883))[2]], [-60, ((diff(x[10](t), t))(0) = HFloat(0.0))[2]], [-58, ((x[11](t))(0) = HFloat(1.71190172933128))[2]], [-56, ((diff(x[11](t), t))(0) = HFloat(0.0))[2]], [-54, ((x[12](t))(0) = HFloat(1.76811705131164))[2]], [-52, ((diff(x[12](t), t))(0) = HFloat(0.0))[2]], [-50, ((x[13](t))(0) = HFloat(1.79644811117089))[2]], [-48, ((diff(x[13](t), t))(0) = HFloat(0.0))[2]], [-46, ((x[14](t))(0) = HFloat(1.79644811117089))[2]], [-44, ((diff(x[14](t), t))(0) = HFloat(0.0))[2]], [-42, ((x[15](t))(0) = HFloat(1.76811705131164))[2]], [-40, ((diff(x[15](t), t))(0) = HFloat(0.0))[2]], [-38, ((x[16](t))(0) = HFloat(1.71190172933128))[2]], [-36, ((diff(x[16](t), t))(0) = HFloat(0.0))[2]], [-34, ((x[17](t))(0) = HFloat(1.62868869443883))[2]], [-32, ((diff(x[17](t), t))(0) = HFloat(0.0))[2]], [-30, ((x[18](t))(0) = HFloat(1.51979026590362))[2]], [-28, ((diff(x[18](t), t))(0) = HFloat(0.0))[2]], [-26, ((x[19](t))(0) = HFloat(1.38692383699642))[2]], [-24, ((diff(x[19](t), t))(0) = HFloat(0.0))[2]], [-22, ((x[20](t))(0) = HFloat(1.23218479067164))[2]], [-20, ((diff(x[20](t), t))(0) = HFloat(0.0))[2]], [-18, ((x[21](t))(0) = HFloat(1.05801345412645))[2]], [-16, ((diff(x[21](t), t))(0) = HFloat(0.0))[2]], [-14, ((x[22](t))(0) = HFloat(0.867156613383085))[2]], [-12, ((diff(x[22](t), t))(0) = HFloat(0.0))[2]], [-10, ((x[23](t))(0) = HFloat(0.66262419483242))[2]], [-8, ((diff(x[23](t), t))(0) = HFloat(0.0))[2]], [-6, ((x[24](t))(0) = HFloat(0.447641796896737))[2]], [-4, ((diff(x[24](t), t))(0) = HFloat(0.0))[2]], [-2, ((x[25](t))(0) = HFloat(0.225599820415749))[2]], [0, ((diff(x[25](t), t))(0) = HFloat(0.0))[2]], [2, ((x[26](t))(0) = HFloat(0.0))[2]], [4, ((diff(x[26](t), t))(0) = HFloat(0.0))[2]], [6, ((x[27](t))(0) = HFloat(-0.225599820415749))[2]], [8, ((diff(x[27](t), t))(0) = HFloat(0.0))[2]], [10, ((x[28](t))(0) = HFloat(-0.447641796896737))[2]], [12, ((diff(x[28](t), t))(0) = HFloat(0.0))[2]], [14, ((x[29](t))(0) = HFloat(-0.66262419483242))[2]], [16, ((diff(x[29](t), t))(0) = HFloat(0.0))[2]], [18, ((x[30](t))(0) = HFloat(-0.867156613383085))[2]], [20, ((diff(x[30](t), t))(0) = HFloat(0.0))[2]], [22, ((x[31](t))(0) = HFloat(-1.05801345412645))[2]], [24, ((diff(x[31](t), t))(0) = HFloat(0.0))[2]], [26, ((x[32](t))(0) = HFloat(-1.23218479067164))[2]], [28, ((diff(x[32](t), t))(0) = HFloat(0.0))[2]], [30, ((x[33](t))(0) = HFloat(-1.38692383699642))[2]], [32, ((diff(x[33](t), t))(0) = HFloat(0.0))[2]], [34, ((x[34](t))(0) = HFloat(-1.51979026590362))[2]], [36, ((diff(x[34](t), t))(0) = HFloat(0.0))[2]], [38, ((x[35](t))(0) = HFloat(-1.62868869443883))[2]], [40, ((diff(x[35](t), t))(0) = HFloat(0.0))[2]], [42, ((x[36](t))(0) = HFloat(-1.71190172933128))[2]], [44, ((diff(x[36](t), t))(0) = HFloat(0.0))[2]], [46, ((x[37](t))(0) = HFloat(-1.76811705131164))[2]], [48, ((diff(x[37](t), t))(0) = HFloat(0.0))[2]], [50, ((x[38](t))(0) = HFloat(-1.79644811117089))[2]], [52, ((diff(x[38](t), t))(0) = HFloat(0.0))[2]], [54, ((x[39](t))(0) = HFloat(-1.79644811117089))[2]], [56, ((diff(x[39](t), t))(0) = HFloat(0.0))[2]], [58, ((x[40](t))(0) = HFloat(-1.76811705131164))[2]], [60, ((diff(x[40](t), t))(0) = HFloat(0.0))[2]], [62, ((x[41](t))(0) = HFloat(-1.71190172933128))[2]], [64, ((diff(x[41](t), t))(0) = HFloat(0.0))[2]], [66, ((x[42](t))(0) = HFloat(-1.62868869443883))[2]], [68, ((diff(x[42](t), t))(0) = HFloat(0.0))[2]], [70, ((x[43](t))(0) = HFloat(-1.51979026590362))[2]], [72, ((diff(x[43](t), t))(0) = HFloat(0.0))[2]], [74, ((x[44](t))(0) = HFloat(-1.38692383699642))[2]], [76, ((diff(x[44](t), t))(0) = HFloat(0.0))[2]], [78, ((x[45](t))(0) = HFloat(-1.23218479067164))[2]], [80, ((diff(x[45](t), t))(0) = HFloat(0.0))[2]], [82, ((x[46](t))(0) = HFloat(-1.05801345412645))[2]], [84, ((diff(x[46](t), t))(0) = HFloat(0.0))[2]], [86, ((x[47](t))(0) = HFloat(-0.867156613383085))[2]], [88, ((diff(x[47](t), t))(0) = HFloat(0.0))[2]], [90, ((x[48](t))(0) = HFloat(-0.66262419483242))[2]], [92, ((diff(x[48](t), t))(0) = HFloat(0.0))[2]], [94, ((x[49](t))(0) = HFloat(-0.447641796896737))[2]], [96, ((diff(x[49](t), t))(0) = HFloat(0.0))[2]], [98, ((x[50](t))(0) = HFloat(-0.225599820415749))[2]], [100, ((diff(x[50](t), t))(0) = HFloat(0.0))[2]]]

 (4)
> 

pointPlot := pointplot(dataPoints, symbol = solidcircle, symbolsize = 8, color = blue)

Error, (in plots:-pointplot) points cannot be converted to floating-point values

  
 

NULL

Download dsolve_fix_3.mw

@Hullzie16 These plots are created ...

Paras31 130
September 11 2024
0 0

@Hullzie16 These plots are created by Matlab for my research. You can read this post and maybe this will help you to understand what I want to create Numerical Simulation of the Discrete Klein-Gordon Equation: A Study on Damped, Driven Nonlinear Wave Systems with Spatially Extended Initial Conditions. 

Thank you very much for your assistance...

Paras31 130
September 10 2024
0 0

@Hullzie16 I tried some ideas and Maple created some useful results, but still, I can not figure out how to create the above results. Thank you very much for your effort!

> 

restart;
with(plots):
> 

L := 200;
K := 99;
kappa := 1;
omegaD := 1;
beta := 1;
delta := 0.05;
j := 2;
tmax := 3000;
h := L/(K + 1);
nsp := [(-L/2 + h*i) $ (i = 0 .. L/h)]:
km := nops(nsp);
omegaD2 := h^2*omegaD^2;
deltaHat := h*delta;
a := 1.8;
var := seq(x[i](t), i = 1 .. km):
initialPositions := seq(x[i](0) = a*sin(j*h*Pi*nsp[i]/L), i = 1 .. km):
initialVelocities := seq(D(x[i])(0) = 0, i = 1 .. km):

200

  

99

  

1

  

1

  

1

  

0.5e-1

  

2

  

3000

  

2

  

101

  

4

  

.10

  

1.8

 (1)
> 

boundaryConditions := [x[1](t) = 0, x[km](t) = 0]

[x[1](t) = 0, x[101](t) = 0]

 (2)
> 

equations := seq(diff(x[n](t), t, t) - kappa*(x[n + 1](t) - 2*x[n](t) + x[n - 1](t)) + deltaHat*diff(x[n](t), t) - omegaD2^2*(x[n](t) - beta*x[n](t)^3) = 0, n = 2 .. km - 1), diff(x[1](t), t, t) = 0, diff(x[101](t), t, t) = 0:
> 

sol := dsolve({equations, initialPositions, initialVelocities}, maxfun=0,numeric, method = rkf45,output=listprocedure)

`[Length of output exceeds limit of 1000000]`

 (3)
> 

# Use odeplot to plot the trajectory of x[50](t)
odeplot(sol, [[t, x[67](t)]], t = 0 .. 100, numpoints = 3500, title = "Trajectory of x[100](t)", axes = boxed, gridlines=true, color = ["Blue"], size = [500, 300]);

 
> 

# Plot x[67](t) vs u[67](t) = diff(x[67](t), t)
odeplot(sol, [[x[67](t), diff(x[67](t), t)]], t = 0 .. 100, numpoints = 3500, title = "x_n vs u_n", axes = boxed, gridlines=true, color = ["Blue"], size = [500, 300]);

 
 

 

Download dsolve_fix_3.mw

Oh ok!! @acer I will have it o...

Paras31 130
September 10 2024
0 0

Oh ok!! @acer I will have it on my mind. I will not do that again

1 2 3 4 Page 1 of 4