Maple 2024 Questions and Posts

These are Posts and Questions associated with the product, Maple 2024

Bernoulli first order ode has form as show in wikipedia  and also on Maple own site as

Notice that it is P(x)*y above and not P(x)* y^(-1) so the y(x) must be linear in that term.   But when I give Maple this ode

ode:=diff(y(x),x) + x*y(x)^(-1)= y(x)^(-1);

Which is clearly not of the form above, it solves it as Bernoulli.  In the above ode, P(x) is x and Q(x) is 1 and n is -1.

The ode advisor correctly said it is separable. But trace shows it used Bernoulli. Also when asking it to solve it as Bernoulli, it does.

What Am I missing here?  Is it not wrong for Maple to use Bernoulli method on this ode which is not Bernoulli?

Worksheet below

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1805 and is the same as the version installed in this computer, created 2024, September 3, 11:35 hours Pacific Time.`

libname;

"C:\Users\Owner\maple\toolbox\2024\Physics Updates\lib", "C:\Program Files\Maple 2024\lib"

restart;

ode:=diff(y(x),x) + x*y(x)^(-1)= y(x)^(-1);
IC:=y(1) = 0;
DEtools:-odeadvisor(ode);

diff(y(x), x)+x/y(x) = 1/y(x)

y(1) = 0

[_separable]

infolevel[dsolve]:=5;

5

dsolve(ode,y(x));  #why this says it solved it as Bernoulli ?

Methods for first order ODEs:

--- Trying classification methods ---

trying a quadrature

trying 1st order linear

trying Bernoulli

<- Bernoulli successful

y(x) = (-x^2+c__1+2*x)^(1/2), y(x) = -(-x^2+c__1+2*x)^(1/2)

dsolve(ode,y(x),[Bernoulli])

Classification methods on request

Methods to be used are: [Bernoulli]

----------------------------

* Tackling ODE using method: Bernoulli

--- Trying classification methods ---

trying Bernoulli

<- Bernoulli successful

y(x) = (-x^2+c__1+2*x)^(1/2), y(x) = -(-x^2+c__1+2*x)^(1/2)

 

 

Download why_this_ode_bernullli_sept_15_2024.mw

Dear Maple community,

I was going to export my math assignment about vector functions as an PDF, when I noticed that it butchered every column vector printed out in blue in the document as math-output. I only have this issue when exporting my work.

How do I make Maple display both coordinates of vector in the blue output field?

I produced the below image as an illustration of my problem:

Thank you in advance.

I am now running eveything in cmaple.exe for number of reasons. But I find the plot generated do not show grid lines as it does in worksheet (GUI). 

Below is worksheet I used, and below that the .mpl file used for command line with the command used.
 

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

restart;

plot(x, x=0..1, axis=[gridlines=[10, color=blue]]);

 


 

Download why_no_grid_in_cmaple_sept_15_2024.mw

Here is the T.mpl file

currentdir("C:/tmp");  #change as needed

#the following commnad below export a plot to a file

plotsetup(ps, plotoutput="t.ps",plotoptions=`color,noborder`);

p0:=plot(x, x=0..1, axis=[gridlines=[10, color=blue]]):

print(p0); #this will create t.ps in same folder

plotsetup(default); #rest back to default

quit;

 

Now I run the above as follows

"C:\Program Files\Maple 2024\bin.X86_64_WINDOWS\cmaple.exe"   T.mpl

Next I convert t.ps to t.pdf using Adobt PDF pro. The output is this

You see, there are no grid lines. I tried changing the color to BLUE, but it made no difference. Tried difference colors also. nothing worked.

Is there a way to keep the grid lines when using cmaple.exe?

I also attached  t.pdf converted by adobe from the .ps file. Mapleprime will not let me attached t.ps file.

t.pdf

 

 

Is this new error in Maple? I do not have older versions to check.

Solved this ode as series. When calling odetest on it, it gives this error as shown in worksheet below.

Maple 2024.1.

Any work arounds?
 

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1805 and is the same as the version installed in this computer, created 2024, September 3, 11:35 hours Pacific Time.`

libname;

"C:\Users\Owner\maple\toolbox\2024\Physics Updates\lib", "C:\Program Files\Maple 2024\lib"

restart;

ode:=-(A/( (1-x)*x^4) - B/( x^4*(1-x)^2)+C/( (1-x)*x^2))*y(x)+1/(1-x)*diff(y(x),x)+diff(y(x),x$2)=0;

-(A/((1-x)*x^4)-B/(x^4*(1-x)^2)+C/((1-x)*x^2))*y(x)+(diff(y(x), x))/(1-x)+diff(diff(y(x), x), x) = 0

Order:=3;
maple_sol:=dsolve(ode,y(x),'series','point'=1);

3

y(x) = c__1*(x-1)^(1-(1-B)^(1/2))*(series(1-((-4*B+A+C)/((1-(1-B)^(1/2))^2+B-1))*(x-1)+((4*A*(1-(1-B)^(1/2))^2-10*B*(1-(1-B)^(1/2))^2+2*C*(1-(1-B)^(1/2))^2+A^2-4*B*A+2*A*C+6*B^2-6*B*C+C^2-4*A+10*B-2*C)/(((1-(1-B)^(1/2))^2+B-1)*((1-(1-B)^(1/2))^2+B+2-2*(1-B)^(1/2))))*(x-1)^2+O((x-1)^3),x = 1,3))+c__2*(x-1)^(1+(1-B)^(1/2))*(series(1-((-4*B+A+C)/((1+(1-B)^(1/2))^2+B-1))*(x-1)+((4*A*(1+(1-B)^(1/2))^2-10*B*(1+(1-B)^(1/2))^2+2*C*(1+(1-B)^(1/2))^2+A^2-4*B*A+2*A*C+6*B^2-6*B*C+C^2-4*A+10*B-2*C)/(((1+(1-B)^(1/2))^2+B-1)*((1+(1-B)^(1/2))^2+B+2+2*(1-B)^(1/2))))*(x-1)^2+O((x-1)^3),x = 1,3))

odetest(maple_sol,ode,'series','point'=1)

Error, (in odetest/series) when calling 'MultiSeries:-multiseries'. Received: 'unable to sort exponents, {(1-B)^(1/2), -(1-B)^(1/2), 1-(1-B)^(1/2), 1+(1-B)^(1/2), 2-(1-B)^(1/2), 2+(1-B)^(1/2), 3-(1-B)^(1/2), 3+(1-B)^(1/2), 4-(1-B)^(1/2), 4+(1-B)^(1/2)}, MultiSeries:-multiseries'

 


 

Download error_odetest_series_sept_14_2024.mw

A classic result states that the equation x3px2qxr=0 with real coefficients p, q, r has positive roots iff p<0, q>0, r<0 and -27r2 - 2p(2p2 - 9q)r + q2(p2 - 4q) ⩾ 0 (see for example this question). 
However, Maple appears unable to find the condition: 

a, b, c := allvalues(RootOf(x^3 + p*x^2 + q*x + r, x), 'implicit'):
RealDomain:-solve({a, b, c} >~ 0, [p, q, r]);
 = 
Warning, solutions may have been lost
                               []

Is there a way to get the above conditions in Maple with as little human intervention as possible (I mean, without a priori knowledge of the theory of polynomials)? 

Edit. An interesting problem is when these three positive roots can further be the lengths of sides of a triangle. For reference, here are some (unenlightening) results from some other software: 

I have number of worksheets open. I set each to use its own engine.

If I have one worksheet running a long loop that takes hrs to complete, or even if I have cmaple.exe running long script in command line, and then at same time I open new worksheet and type

             Physics:-Version(latest)

Then I always get this error

This is described here

It seems because Maple was busy running another worksheet. I never see this error if I am not running any other thing in Maple at the time. ( I actually try not to update Physics while running something else in Maple, but sometimes I do not notice that I have something else running).

When I wait for my other worksheet to finish, then close all of Maple and reopen Maple, I see the latest version is installed, at least this is what Physics:-Version() now says.

My question is:  Is it safe to do Physics:-Version(latest) even when Maple is busy running  computation in other separate worksheet or even running cmaple.exe from windows command line? Is there any concurrency/locking issue on the Physics library being updated while it is being open by another worksheet or another server.exe running in background?

All this is on windows 10 pro.

I'm trying to solve a system of coupled differential equations numerically, but I'm getting the following error

Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations

The error occurs at the dsolve step, despite trying to ensure that all equations and conditions are in the correct form (sets/lists).

Could someone help me identify what I'm missing here?

Thanks in advance!

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

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

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)

sol := dsolve({equations, boundaryConditions, initialPositions, initialVelocities}, var, numeric, method = rkf45, range = 0 .. tmax)

Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations

 

NULL

Download dsolve_error.mw

Why Maple returns y(x)=0 as solution to this ode when even odetest do not validate the solution? Solution satisfies the ode but not the given BC.

Is this a known bug?

I've updated the worksheet now to include few more examples, all using same ode but with different BC's. In all 4 examples, Maple gives solutions when there should not be solution and in one case gives solution which can not even be used for odetest. These are random BC's generated just for testing.
 

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1804. The version installed in this computer is 1802 created 2024, September 3, 11:35 hours Pacific Time, found in the directory C:\Users\Owner\maple\toolbox\2024\Physics Updates\lib\`

libname;

"C:\Users\Owner\maple\toolbox\2024\Physics Updates\lib", "C:\Program Files\Maple 2024\lib"

restart;

#EXAMPLE 1

ode:=diff(diff(y(x), x), x) + 2*diff(y(x), x)/x - y(x) = 0;
IC:=y(-1)=0,D(y)(-infinity)=-1/exp(1);

diff(diff(y(x), x), x)+2*(diff(y(x), x))/x-y(x) = 0

y(-1) = 0, (D(y))(-infinity) = -1/exp(1)

maple_sol:=dsolve([ode,IC]);

y(x) = 0

odetest(maple_sol,[ode,IC])

[0, 0, -exp(-1)]

restart;

#EXAMPLE 2

ode:=diff(diff(y(x), x), x) + 2*diff(y(x), x)/x - y(x) = 0;
IC:=y(infinity)=0,D(y)(-infinity)=-1/exp(1);

diff(diff(y(x), x), x)+2*(diff(y(x), x))/x-y(x) = 0

y(infinity) = 0, (D(y))(-infinity) = -1/exp(1)

maple_sol:=dsolve([ode,IC]);

y(x) = 0

odetest(maple_sol,[ode,IC])

[0, 0, -exp(-1)]

restart

#EXAMPLE 3

ode:=diff(diff(y(x), x), x) + 2*diff(y(x), x)/x - y(x) = 0;
IC:=y(-infinity)=0,D(y)(-infinity)=-1/exp(1);

diff(diff(y(x), x), x)+2*(diff(y(x), x))/x-y(x) = 0

y(-infinity) = 0, (D(y))(-infinity) = -1/exp(1)

maple_sol:=dsolve([ode,IC]);

y(x) = signum((sinh(x)+cosh(x))/x)*infinity

odetest(maple_sol,[ode,IC])

Error, (in signum) signum is not differentiable at 0

#EXAMPLE  4

restart;

ode:=diff(diff(y(x), x), x) + 2*diff(y(x), x)/x - y(x) = 0;
IC:=y(-infinity)=0,D(y)(infinity)=-1/exp(1);

diff(diff(y(x), x), x)+2*(diff(y(x), x))/x-y(x) = 0

y(-infinity) = 0, (D(y))(infinity) = -1/exp(1)

maple_sol:=dsolve([ode,IC]);

y(x) = 0

odetest(maple_sol,[ode,IC])

[0, 0, -exp(-1)]

 


 

Download wrong_sol_dsolve_sept_10_2024.mw

Update sept 13, 2024

Here is one more example I found. In this ode, the IC given leads to division by zero in the solution. Yet Maple for some reason removes the offending part of the solution with that constant, and returns the rest. 

Is one allowed to do this? odetest also do not validate the solution. There should not be solution returned in this case, since it is not possible to find values for the constants of integration given these initial conditions. Right?

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

restart;

ode:=x^2*diff(y(x),x$2)+x*diff(y(x),x)-9*y(x)=0;
IC:=y(0)=1,D(y)(0)=0;
maple_sol_no_IC:=dsolve(ode);

x^2*(diff(diff(y(x), x), x))+x*(diff(y(x), x))-9*y(x) = 0

y(0) = 1, (D(y))(0) = 0

y(x) = c__1*x^3+c__2/x^3

maple_sol_with_IC:=dsolve([ode,IC]); #maple just removed the c2/x^3 part. Since at x=0 gives problem?

y(x) = c__1*x^3

odetest(maple_sol_with_IC,[ode,IC]); #shows solution does not verify

[0, 1, 0]

 


 

Download another_example_solution_given_when_none_exist.mw

 

 

I wanted to change  eq:= 1/2 * sqrt(-2*lambda)  to 1/2 %* sqrt(-2*lambda)  using a rule.

It works outside of rule ofcourse. But when I put %* in the RHS of the rule, maple hangs. It seems it is going into infinite loop.

I tried the trick of using '%*' but this gives syntax error.

Same problem happens when using %. and not just %*

Is there a workaround?

Attached worksheet. Make sure to save all work before trying it.
 

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

Physics:-Version()

`The "Physics Updates" version in the MapleCloud is 1804. The version installed in this computer is 1802 created 2024, September 3, 11:35 hours Pacific Time, found in the directory C:\Users\Owner\maple\toolbox\2024\Physics Updates\lib\`

libname;

"C:\Users\Owner\maple\toolbox\2024\Physics Updates\lib", "C:\Program Files\Maple 2024\lib"

restart;

eq:= 1/2 * sqrt(-2*lambda)

(1/2)*(-2*lambda)^(1/2)

eq:= 1/2 %. sqrt(-2*lambda); #no problem

`%.`(1/2, (-2*lambda)^(1/2))

eq:= 1/2 %* sqrt(-2*lambda); #no problem

`%*`(1/2, (-2*lambda)^(1/2))

restart;

eq:= 1/2 * sqrt(-2*lambda)

(1/2)*(-2*lambda)^(1/2)

applyrule(sqrt(x::anything)/y::anything = 1/y %. sqrt(x),eq); #why this hangs?

 

restart;

eq:= 1/2 * sqrt(-2*lambda)

(1/2)*(-2*lambda)^(1/2)

applyrule(sqrt(x::anything)/y::anything = 1/y %* sqrt(x),eq); #why this hangs?

 


 

Download maples_hangs_applyrule_sept_10_2024.mw

Hello!

I present a simple work-up of rolling a plane curve along a fixed plane curve in 2d space. Maple sources are attached.

Kind regards!

Source.zip

Hello,

Please help me solve this type of equations, I have more of them, but there is only one in the attached file. fsolve does not work. Please help.

Best regards!

tg_Test.mw

The help page of interface('longdelim') states: 

If true, Maple control structures such as if, do, proc, and so on are displayed with the newer-style long ending delimiters such as end if, end do, end proc, and so on. If false, ending delimiters are displayed as fi, od, end, and so on

If I set interface('longdelim' = false):, the Maple Input

try f:=0 catch:end try;

can be converted into 

via , but why is

use f=0 in f+1 end use;

still converted into

 instead of “use f = 0 in f + 1 end”? 

With the attached task I would like to learn how Maple handles polynomials and plots. For this I have chosen a task with an interesting history from 1593. Adriaan van Roomen posed it to F. Vieta, who solved it after a short consideration.

My questions are:
1.) How are very long terms entered and displayed in Maple in an appropriate manner?
2.) How are the graphs of several functions displayed in the same coordinate system?
3.) Which graphics settings must be selected for differentiable functions in order to obtain nice, rounded curves rather than angular ones from the numerical result?
4.) Can the apparently random oscillations of the curve at the end of the interval be suppressed?
5.) The curves for p(x) and sinnx(x) are theoretically identical, since p(x) is the trigonometric expansion of sinnx(x). Their graphs must therefore be identical. How can this be displayed?
6.) The polynomial has many real zeros. How can the zeros be clearly presented in a table?

AF_20240909.mw

This is an excersise in one of Mathematica's tutorials. The problem is to find from all persumtations of [1,2,3,4] those lists which has 2 in them before 3.

But the idea is that 2 can be anywhere before the 3, with possibly zero or more numbers between.

Will show how to do this in that other software, and ask if there is a way using Maple pattern matching or better way to do this in Maple better than what I did.

L=Permutations[{1,2,3,4}]
Cases[L,{___,2,___,3,___}]

gives

{{1, 2, 3, 4}, {1, 2, 4, 3}, {1, 4, 2, 3}, {2, 1, 3, 4}, 
{2, 1, 4,  3}, {2, 3, 1, 4}, {2, 3, 4, 1}, {2, 4, 1, 3}, 
{2, 4, 3, 1}, {4, 1,  2, 3}, {4, 2, 1, 3}, {4, 2, 3, 1}}

We see in all of the above, 2 is before 3.

In that other software, the ___ means there is zero or more things.  I could not find how to do this in Maple's pattern matching. So had to use has and then find the index of 2 and 3 in each list and check if the index of 2 is smaller than the index of 3. Which is kinda awakard and not as elegent as using a pattern, but it gives same result.

L:=combinat:-permute([1,2,3,4]);
map(X->`if`(has(X,2) and has(X,3) and ListTools:-Search(2,X)<ListTools:-Search(3,X),X,NULL) ,L);

Gives

[[1, 2, 3, 4], [1, 2, 4, 3], [1, 4, 2, 3], [2, 1, 3, 4], 
[2, 1, 4, 3], [2, 3, 1, 4], [2, 3, 4, 1], [2, 4, 1, 3], 
[2, 4, 3, 1], [4, 1, 2, 3], [4, 2, 1, 3], [4, 2, 3, 1]]

Can you suggest a way using either patmatch or applyrules in Maple to do the same?  

Could this be done easier than what I did using evalindents without the need for pattern?

I find patterns easier to work with myself.

ps. converting each list to string, and then using regular expression or string matching, is not what I am looking for. 

ps. the check in Maple code of has(X,2) and has(X,3) is not really needed here, since we know each list will have 2 and 3 in them. But I kept these to make it more general for other type of problems where these extra checks might be needed.

I am trying to create the plot for the dynamics of the discrete Klein - the equation gives Gordon equation (DKG) with friction for t=0 and t=3000.

As you can see below, creating the plot for t=0 is pretty easy. However, for t=3000, I tried to create a procedure to produce the desired plots, but when I ran the code, nothing happened. The worksheets are evaluated all the time. There are no results.

I am studying  for my thesis  the dynamics of the discrete Klein - the equation gives Gordon equation (DKG) with friction:

 

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 potential 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 is the discretization parameter. We denote by the distance between the oscillators of the lattice. The chain (DKG) contains linear damping with a damping coefficient, while β s the coefficient of the nonlinear cubic term.

 

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__Error loading this structure.`-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

 

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) "

 

• 

t=0

 

restart; with(plots); a := 2; j := 2; L := 200; K := 99; h := L/(K+1); n := Vector([`$`(-(1/2)*L+h*i, i = 0 .. K+1)]); U_n0 := a*map(proc (x) options operator, arrow; sin(j*Pi*x/L) end proc, h*n)

plot([seq([n[i], U_n0[i]], i = 1 .. K+2)], style = point, color = blue, axes = boxed, gridlines, grid = [true, true], labelfont = [TIMES, 12])

 

Now I am trying to create plot for different values of a.

a = {2, 1.82, 1.85, 1.9, 1.95}, K = 99, L = 200, beta = 1, delta = 0.5e-1, `&omega;__d`^2 = 1

restart; with(plots)

L := 200; K := 99; j := 2; omega_d := 1; beta := 1; delta := 0.5e-1

h := L/(K+1); n := Vector([`$`(-(1/2)*L+h*i, i = 0 .. K+1)])

dt := 0.5e-1; tmax := 3000; tspan := [seq(0+dt*i, i = 0 .. trunc(tmax/dt))]

NULL

a := 1; U0 := a*map(proc (x) options operator, arrow; sin(j*Pi*h*x/L) end proc, n); U0[1] := 0; U0[K+2] := 0; Udot0 := Vector(K+2, 0)

discrete_laplacian := proc (U) local i, Laplacian; Laplacian := Vector(K+2); for i from 2 to K+1 do Laplacian[i] := U[i+1]-2*U[i]+U[i-1] end do; Laplacian[1] := 0; Laplacian[K+2] := 0; return Laplacian end proc

NULL

U := U0; Udot := Udot0; for t to num_steps do Lap := discrete_laplacian(U); Ucube := map(proc (x) options operator, arrow; x^3 end proc, U); for i from 2 to K+1 do Udot[i] := Udot[i]+dt*(2*omega_d^2*U[i]-beta*Ucube[i]-delta*Udot[i]+Lap[i]) end do; for i from 2 to K+1 do U[i] := dt*Udot[i]+U[i] end do; U[1] := 0; U[K+2] := 0; if `mod`(t, 1000) = 0 then plotU := plot(n, U, style = line, color = blue, title = cat("Displacement at t=", t*dt), labels = ["n", "Displacement"], grid = [true, true]); display(plotU) end if end do

NULL


My target is the following plots

Download Numerical_Simulation_of_discrete_Klein-Gordon_Equation.mw

1 2 3 4 5 6 7 Last Page 1 of 19