Paras31

140 Reputation

3 Badges

0 years, 202 days
Hellenic Open University
Mathematician

Social Networks and Content at Maplesoft.com

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

@mmcdara

I intended to use delta=0 for a general case and compare it with the result of the first worksheet. I wanted to see if they would have the same behavior the following

ODE__1 := m*diff(x(t), t, t) + k*x(t) = 0

when I set all other parameters equal to zero in the following equation

ODE__2 := m*diff(x(t), t, t) + delta*diff(x(t), t) + k*x(t) = 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 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


 

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

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

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

@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

1 2 3 4 Page 1 of 4