## 13598 Reputation

19 years, 197 days

## Pendulum...

I notice that the equation eqt considered as giving x implicitly as a function of t solves the well-studied pendulum equation:

diff(x(t), t, t) + k*sin(x(t)) =0;

For information you could go to
https://en.wikipedia.org/wiki/Pendulum_(mechanics)

So do the solving numerically:

restart;
ode := diff(x(t), t, t)+ k*sin(x(t));
res:=dsolve({ode,x(0)=a,D(x)(0)=0},numeric,parameters=[a,k]);
res(parameters=[1,.3]); #Setting parameters
plots:-odeplot(res,[t,x(t)],-20..20);
res(parameters=[1,1.5]); #Setting new parameters
plots:-odeplot(res,[t,x(t)],-20..20);


Be careful: If you reexecute the first plotting command then the most recent parameters are used.
You can, however, save the plots as in this version:

restart;
ode := diff(x(t), t, t)+ k*sin(x(t));
res:=dsolve({ode,x(0)=a,D(x)(0)=0},numeric,parameters=[a,k]);
res(parameters=[1,.3]); #Setting parameters
plots:-odeplot(res,[t,x(t)],-20..20,thickness=3); p1:=%:
res(parameters=[1,1.5]); #Setting new parameters
plots:-odeplot(res,[t,x(t)],-20..20,color=blue); p2:=%:
plots:-display(p1,p2);


## Selecting the input...

Using the file provided I tried selecting the cosine expression. Then right clicked on that and chose convert to 1D-Math input.
That worked. The same worked on sqrt(4).

## A suggestion...

I don't know what you mean by wanting to create a vector x.
But in order to find out you could look at this code:

restart;
n:=3:
Theta:=<0.2,0.3,1.5>;
A:=<2,6,3>;
F:=<10,4,7>;
f(.9);
plot(f(x),x=-3..10);


Guessing: Could this continuation be it?

N:=10000:
X:=Vector(N,i->0.001*i);
Y:=f~(X);
plot(X,Y); # Just for illustration


## pdsolve...

restart;
pde:=diff(<u1(x,t),u2(x, t)>, t) = Matrix([[0, mu*k^2/2], [2*A^2 - mu*k^2/2, 0]]).<u1(x,t),u2(x, t)>;
pdsolve(pde);


## Use subs...

You could use subs as in this code:

A:=Int(r^2, m);
B:=subs(m=(m=a..b),A);


## simplify...

On my Maple 2022.2 I just added simplify:

restart;
g := piecewise((1/2)*t < 1/2, 0, (1/2)*t < 3/2, 1, 0)*piecewise(z/t < 1, 0, z/t < 2, 10080-60480*z/t+156240*z^2/t^2-226800*z^3/t^3+202230*z^4/t^4-113400*z^5/t^5+39060*z^6/t^6-7560*z^7/t^7+630*z^8/t^8, 0)/t;

Int(g, t=1..3);
value(%);
simplify(%);


That brought an explicit answer. I don't have Maple 2015 on the computer I have available here, but it works on Maple 2021.2 too.

## _bigPi...

You can see the code for evalf/constant/Pi by doing showstat(evalf/constant/Pi);

There you see that there are 4 cases.
The global variable called _bigPi initially has 1000 digits.
Try this:

restart;
showstat(evalf/constant/Pi);
length(op(1,_bigPi)); # 1000
_bigPi;
evalf[2000](Pi);
length(op(1,_bigPi)); # 10000
_bigPi;


If you continue the session above with evalf[10001](Pi) you will find that now
length(op(1,_bigPi)) = 10001.
After a restart it will be back at 1000.

## The second loop...

You may have expected something else than what comes out of the following slightly corrected version, if so let us know:

restart;
n:=5;
R:=Matrix(5,n);

for i from 1 to n do
R[..,i]:=sin(i*x);
end do;
R;
for i from 1 to 5 do
#i:=i; # This line doesn't do anything useful, but isn't harmful
x:=i;
R[i,..] := sin(x);
end do;
R;


## Ignore the first four...

dsolve assumes that c and g are constants, so the first 4 solutions are constant solutions and 2 of these are imaginary when c and g are positive, which you assume.

I assume that you aren't interested in constant solutions. Thus the fifth solution is the one you want.

I don't think you can convert the JacobiSN function into a trigonometric function.

## The local variable Omega~...

Here is an attempt to explain the situation when using assume.MaplePrimes22-09-29_assume.mw

 > restart;
 > A := Omega;          # A now evaluates to Omega assume(0 < Omega);   # Omega will now evaluate to the local variable Omega~ B := Omega;          # B will now evaluate to Omega~ by the default full evaluation: first Omega then Omega~

 > type(Omega,local); # true

 > # All of the following 'is' commands result in 'true':
 > is(A>0);

 > is(Omega>0);

 > is(0 < B);

 > is(B=Omega);

 > is(A=B);

 > is(A=Omega);

Originally Omega, renamed Omega~:
is assumed to be: RealRange(Open(0),infinity)

Originally Omega, renamed Omega~:
is assumed to be: RealRange(Open(0),infinity)

Originally Omega, renamed Omega~:
is assumed to be: RealRange(Open(0),infinity)

 > Omega:='Omega';      # This assignment makes Omega just evaluate to its own name Omega

Omega:

 > about(A);            # A will first evaluate to Omega and Omega now just evaluates to its own name

Omega:

 > about(B);            # B was assigned to Omega when full evaluation gave Omega~, thus B evaluates to Omega~.

Omega:
is assumed to be: RealRange(Open(0),infinity)

 > is(B,local);       # true

 >

## unwith and forget...

@acer I notice in Carl's second example that unwith followed by forget will restore the default diff.
That doesn't explain what's going on, but at least a restart is not necessary.

## Use normal...

restart;
expr:=A/B+C;
#combine(expr);
normal(expr);


## gamma...

In Maple gamma is Euler's constant (see ?gamma), which is approximately equal to 0.5772156649. Try evalf(gamma).
Either of the remedies proposed by vv or tomleslie will work, however, in your case since either one is unaffected by that fact as only alpha*gamma appears. I also pointed that out in my answer to your previous question dealing with the same system. (https://mapleprimes.com/questions/234300-How-Can-I-Fix-Error-in-Dsolvenumericprocessinput ).
In general if you want to use gamma as just another name like alpha and beta, you can declare gamma local at the top of your worksheet as in:

restart;
local gamma;
## The rest ......

This only works in recent Maple versions. The one you use, Maple 2018, is included.

## Several typos and other problems...

I deleted all the output that should not have been copied as Tom Leslie says.

Besides I corrected some typos and errors. So in that process I may have made my own errors.
The following appears to be what you want. Notice in particular that alpha*gamma never appears in your original equations, but alpha*m does. I changed that product to a single name: am.

restart;

with(plots):

U := [W(t), X(t), Y(t), Z(t)];
f := (W, X, Y, Z) -> A-mu*W-delta1*Y-delta2*Z;
g := (W, X, Y, Z) -> beta*(W-X-Y-Z)*X-mu*X-X;

h := (W, X, Y, Z) -> am*X-(mu+delta1)*Y;

i := (W, X, Y, Z) ->(-am+1)*X-(mu+delta2)*Z;

sys1 := diff(U, t) =~ [f(op(U)), g(op(U)), h(op(U)), i(op(U))];

ic1 := W(0) = W__0, X(0) = X__0, Y(0) = Y__0, Z(0) = Z__0;

equil := solve({f(w0, x0, y0, z0) = 0, g(w0, x0, y0, z0) = 0, h(w0, x0, y0, z0) = 0, i(w0, x0, y0, z0) = 0}, {w0, x0, y0, z0});

param1 := {A = 1, beta = 0.9e-3, delta1 = 0.139e-2, delta2 = 0.134e-2, mu = 0.132e-2, W__0 = 750, X__0 = .1, Y__0 = .2, Z__0 = .6, am = 0.4e-2};

param2 := {A = 1, beta = 0.9e-3, delta1 = 0.139e-2, delta2 = 0.134e-2, mu = 0.132e-2, W__0 = 755, X__0 = .2, Y__0 = .4, Z__0 = .9, am = 0.4e-2};

param3 := {A = 1, beta = 0.9e-3, delta1 = 0.139e-2, delta2 = 0.134e-2, mu = 0.132e-2, W__0 = 760, X__0 = .3, Y__0 = .6, Z__0 = 1, am = 0.4e-2};

param4 := {A = 1, beta = 0.9e-3, delta1 = 0.139e-2, delta2 = 0.134e-2, mu = 0.132e-2, W__0 = 765, X__0 = .4, Y__0 = .7, Z__0 = 1.4, am = 0.4e-2};

sol1 := dsolve(eval({ic1, op(sys1)}, param1), U, numeric);

odeplot(sol1,[t,X(t)],0..30);
for j to 4 do
p[j]:=odeplot(sol1,[t,U[j]],0..30);
end do:
display(Array([seq(p[j],j=1..4)]));


The plot of X(t):
MaplePrimes22-06-06_odesys.mw

## subs or Physics:-Assume...

This is tricky. That is why Physics:-Assume exists, I believe.

restart;
assume(A < B);
S := 2/(B-A);
A;
subs(A=5,B=10,S); # 2/5
################
restart;
Physics:-Assume(A < B);
S := 2/(B-A);