## 11973 Reputation

16 years, 57 days

## Question (c) and (d)...

```restart;
## The right hand sides of your system:
f:=(x,y)->-x+x^2-3*x*y;
g:=(x,y)->-y-5*x*y+y^2;
## The proposed Lyapunov function:
V := (x,y)->x^2/2 + y^2/2;
## Clearly V(0,0) = 0 and V is positive for (in fact) all (x,y)<>(0,0).
## Now consider:
dV:=<diff(V(x,y),x),diff(V(x,y),y)>.<f(x,y),g(x,y)> assuming real;
## For local stability at (0,0) we need only consider the lowest order terms of dV.
## Those are of order 2:
mtaylor(dV,[x,y],3); # -x^2-y^2: Clearly negative for (x,y)<>(0,0)
```

Conclusion: dV is negative in a neighborhood U of (0,0) if (x,y) <>(0,0).
Thus (0,0) is (locally) asymptotically stable.
To give an explicit neighborhood we can make a few estimates:

```dV2:=mtaylor(dV,[x,y],3);
dV3:=expand(dV-dV2);
dV3a:=map(abs,dV3);
## Crude, but correct:
dV3a<=(abs(x)+3*abs(y)+5*abs(x)+abs(y))*(x^2+y^2);
```

Thus if
6*abs(x)+4*abs(y) <= k and k < 1;
we have dV = dV3+dV2 <= dV3a+dV2 <=
(6*abs(x)+4*abs(y))*(x^2+y^2)+dV2 = -(1-k)*(x^2+y^2) , so dV < 0 if also (x,y)<>(0,0).

So a neighborhood U that works is defined by the inequality: 6*abs(x)+4*abs(y) <= k, where k < 1.

The basin of attraction of (0,0) is bounded by the cyan colored separatrices to the 3 saddle points in this phase plane image:

## x and X and more...

Since you have

`X:=x(tau): `

you cannot use X in your IC. You must use x.
Furthermore w1 needs to be given a value for the numerical solution to work:

```IC:=x(0)=0,D(x)(0)=0;
##
sys1:={sys,IC}:
dsolve(eval(sys1,w1=1),numeric,method=rkf45); # Example
```

Whether a symbolic solution ever appears you will see if you are more patient than I am.

## Check against dsolve/numeric...

You can check your results with dsolve/numeric.
If you remove the method option method=classical[rk4], then you get the default method=rkf45 (then also remove the stepsize option). Using the parameters option allows you to change the values of x0, k, and f in a very simple manner as below.

```## Assuming that m, F, and Omega have been defined as in your worksheet.
k:='k': f:='f':
### Using method=classical[rk4] for direct comparison
solRK4:=dsolve({ode,x(0)=x0,D(x)0)=0},numeric,method=classical[rk4],stepsize=0.001,parameters=[x0,k,f]);
solRK4(parameters=[3,0.1,0]); #Setting parameters
plots:-odeplot(solRK4,[x(t),diff(x(t),t)],0..10);
plots:-odeplot(solRK4,[t,1/2*diff(x(t),t)^2+V(x(t))],0..10); #Energy plot
plots:-odeplot(solRK4,[[t,x(t)],[t,diff(x(t),t)]],0..10,size=[800,default]);
```

You can set parameters as done above as you wish.

## depends...

You may try depends.

`depends(expr,a);`

where expr is the polynomial or equation.

## Singularity...

There could very likely be a singularity, not only a numerical problem.
Digits:=15:
# In the loop used this try..catch structure:

```try
sh0:=SH(0):
sv0:=SV(0):
Lx0:=LA(0):
print('i'=i);
catch:
print("Breaking",'i'=i); break
end try;
```

around the 3 assignments.
Furthermore, in the dsolve command I put in stiff=true, maxfun=10^6, and removed the range argument.
Then the loop ran just once. The break came at i = 2 probably because of an actual singularity strongly indicated here:

```sol(parameters);
odeplot(sol,[x,L(x)],4.238507..H-He);
odeplot(sol,[x,sh(x)],4.238507..H-He);
odeplot(sol,[x,sv(x)],4.238507..H-He);
odeplot(sol,[x,u(x)],4.238507..H-He);
```

## Execute again...

Looking at the output from the procedure definition of S1 (the blue stuff) I see psi inside the procedure body although the input appears to have phi (actually varphi).
My guess is that you had a first version defined like that. Then you realized that that was wrong and changed the input. You didn't execute the revised definition, however.

## Remove 'numeric'...

Just remove 'numeric':

```int(piecewise(0 < -k5*x1*x2+k1*x1 and -k5*x1*x2+k1*x1 < 1, 1, -k5*x1*x2+k1*x1 <= 0 and 1 <= -k5*x1*x2+k1*x1, 0)*piecewise(0 < -k5*x1*x2+k2*x2 and -k5*x1*x2+k2*x2 < 1, 1, -k5*x1*x2+k2*x2 <= 0 and 1 <= -k5*x1*x2+k2*x2, 0), k1 = 0 .. 1, k2 = 0 .. 1, k5 = 0 .. 1, x1 = 0 .. 1, x2 = 0 .. 1);
int(piecewise(0 < x+y and x+y < 1, 1, x+y <= 0 and 1 <= x+y, 0)*piecewise(0 < y*x^2 and y*x^2 < 1, 1, y*x^2 <= 0 and 1 <= y*x^2, 0), x = 0 .. 1, y = 0 .. 1);

That numeric integration takes a long time or has its problems is likely due to your integrands not being continuous.

## Nonsense...

Maybe it is a slight exaggeration to call your problem nonsense.
But my reason is this. Your system can in vector form be written as
X'(t) = f(t, X(t)), X(t0) = X0  (*)
where X and X0 are vectors, in your case of dimension 2, and f is a vector valued function.
The initial value problem (*) is equivalent to the integral equation
X(t) = X0 + int( f(s, X(s) ), s=t0 .. t);
As the value of an integral stays the same if the integrand is changed on a set of measure zero, your problem is equivalent to the simple problem

`sys:={diff(x(t),t)=y(t), diff(y(t),t)=-x(t), x(0)=0,y(0)=1};`

## A missing matching parenthesis...

In the definition of c a final parenthesis ')'  is missing. After that no complaint from the parser.

## Singularities...

Isn't it a problem that your function has infinitely many singularities on the positive real axis
(s[n], n=1..infinity) where s[n] -> infinity as n - > infinity
?
To see that there are singularities as described (and as expected because of the tan function) you could do:

```restart;
F:=1/(s^2*(1+(1-1/(2+2*s))/s)*(1+tan((1/2)*Pi/sqrt((1-1/(2+2*s))/s))/sqrt((1-1/(2+2*s))/s)));
F2:=normal(convert(F,sincos));
A:=op([1,1],indets(F2,specfunc(sin)))/Pi;
limit(A,s=infinity);
sp1,sp2:=solve(A=p,s);
evalf(seq(sp1,p=1..10));
dF2:=denom(F2);
seq(fsolve(dF2=0, s=eval(sp1,p=n)..eval(sp1,p=n+1)),n=0..9);
```

The sine term in the denominator dominates the cosine term when s > 1 and tends to zero as s->infinity. Thus the singularities will be close to those values of s that make the argument to sine an integral multiple of Pi. Those values are given by sp1 for p = 1..infinity.
The first 10 actual singularities are found in the last command by fsolve.

## Your example as an initial value problem...

Thanks to acer for suggesting major parts of the following approach.

Making use of dsolve 20000 times takes a while. So if the system could be turned into an initial value problem we can use the parameters option.

Although your system is really a boundary value problem since it has four functions with conditions given at H-He and one function (u) with u(0) = 0 we can consider instead posing the condition u(H-He)=0.
Your system is special. The ode for u has the form diff(u(x),x) = f( Jir(x) ) = g(x), i.e. it doesn't depend on u(x).
Thus the u(x) is just obtained by integration, u(x) = int(g(xi), xi = H-He .. x) (not that we shall do that literally).
The function u we actually want (ua, say) is then given by ua(x) = u(x)-u(0).
I have modified other parts, but that is really not important.

```restart;
Digits:=15:
phi0:=0.72:
rho0:=646.8:
E:=100*10^6:
nu:=1/3:
pc0:=0.3*10^6:
pv0:=0.6*10^6:
g:=9.80665:
Mp:=rho0*100/(10^6*365*24*60*60):
lambda:=E*nu/((1+nu)*(1-2*nu)):
mu:=E/(2*(1+nu)):
K:=lambda+2/3*mu:
eta:=100*10^6*10^6*365*24*60*60:
a:=1.15:
A:=-2/3+a*sqrt(3)/3:
B:=-1/3-a*sqrt(3)/3:
m:=5:
n:=3:
Ts:=60*10^6*365*24*60*60:
Tf:=2*Ts:
####
Md:=Mp*t:
phi:=1-(1-phi0)/Jir(x);
pc:=pc0*(ln(phi)/ln(phi0))^m;
mp:=-(m*(1-phi0)*pc)/(Jir(x)*phi*ln(phi));
Gp:=-(K+2*sqrt(3)/3*mu*a)/(K+a^2*mu+mp);
beta11:=K-2/3*mu;
beta21:=K+4/3*mu;
beta12:=(K-sqrt(3)/3*mu*a)*Gp+beta11;
beta22:=(K+2*sqrt(3)/3*mu*a)*Gp+beta21;
####
Delta1:=1/(rho0*g);
H:=beta21*Delta1*(1-exp(-g/beta21*Md));
Lambda:=1-rho0*g/beta21*(H-x);
####
sh:=beta11*ln(Lambda);
sv:=beta21*ln(Lambda);
Le:=evalf(exp(-pc0/(K+2*sqrt(3)/3*mu*a)));
Jire:=1;
He:=beta21*Delta1*(1-Le);
Te:=fsolve(H=He,t);
she:=eval(sh,{x=0,t=Te});
sve:=eval(sv,{x=0,t=Te});
ue:=-Mp*g/beta21*He;
unassign('Lambda','sh','sv','H'):
###############################
sys_ode:=diff(sv(x),x)=rho0*g/Lambda(x),diff(sh(x),x)=beta12/beta22*rho0*g/Lambda(x),diff(Lambda(x),x)=rho0*g/beta22,diff(u(x),x)=-Mp*g/beta22,diff(Jir(x),x)/Jir(x)=-Gp*diff(Lambda(x),x)/Lambda(x);
### Your ics has u(0)=0. The new ones:
icsNew:={sv(H-He)=sve,sh(H-He)=she,Lambda(H-He)=Le,Jir(H-He)=Jire,u(H-He)=0};
resNew:=dsolve({sys_ode} union icsNew,numeric,[sv(x),sh(x),Lambda(x),Jir(x),u(x)],parameters=[H],
output=listprocedure,abserr=1e-14,relerr=1e-12); #optimize may help a little
SV,SH,JIR,U:=op(subs(resNew,[sv(x),sh(x),Jir(x),u(x)]));
resNew(parameters=[He]); #Setting the parameter H to He for a test:
SV(0),SH(0),JIR(0),U(0);
#######
nn:=20000:
CC1:=Matrix(datatype=float[8]):
CC2:=Matrix(datatype=float[8]):
TR1:=Matrix(datatype=float[8]):
TR2:=Matrix(datatype=float[8]):
EN1:=Matrix(datatype=float[8]):
EN2:=Matrix(datatype=float[8]):
####
t0:=time():
dt:=(Ts-Te)/nn:
t:=Te:
HH:=He:
for i from 1 to nn+1 do
resNew(parameters=[HH]); #Setting the parameter H to HH
CC2(i,1..2):=Array([t/(10^6*365*24*60*60),HH]):
sh0:=SH(0):
sv0:=SV(0):
TR2(i,1..2):=Array([-1/3*sh0-2/3*sv0,evalf(sqrt(3)/3*(sh0-sv0))]):
Jir0:=JIR(0):
phix0:=1-(1-phi0)/Jir0:
EN2(i,1):=pc0*(ln(phix0)/ln(phi0))^m:
up:=U(HH-He)-U(0); # Since U is just an integral.
uh:=ue+up:
fvp:=evalf[30](A*sh0+B*sv0-(pv0*(ln(phix0)/ln(phi0))^n)):
if fvp>=0 then
break
else
dH:=dt*(Mp/rho0+uh):
HH:=HH+dH:
t:=t+dt:
end if:
end do:

time()-t0;
i;
```

The time used on my computer was 3.140s and i evaluates to 742.
MaplePrimes18-08-21_odesys_inits.mw

## _MonteCarlo...

Using method= _MonteCarlo and given a large enough epsilon:

`int(add(cat(x,i),i=1..10)^2,seq(cat(x,i)=0..1,i=1..10)],numeric,method=_MonteCarlo,epsilon=0.5e-4);`

I got 25.83247877.

## A dsolve/numeric workaround...

First an observation: In Maple 8 the 3 results of the evalf command in the following have no imaginary parts at all:

```M:=MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi);
seq(evalf(eval(M,x=xx)),xx=[0.4,0.5,0.6]);
```

In Maple 12, 15, 16, ..., 2018 there are very small imaginary parts in all 3.
This explains the difference in the results of this workaround for the integral:

```u:=(-5*ln(x)^4*Pi^4-20*ln(x)^2*Pi^4-8*Pi^4+120*MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi))/(120*Pi^4*(-1+x)^2);
ode:=diff(y(x),x)=u;
res:=dsolve({ode,y(4/10)=0},numeric); #OK in Maple 8, but not in 12, 15...2018
##A workaround for the workaround in 12, 15...2018
res:=dsolve({ode,y(4/10)=1e-22*I},numeric);
res(6/10);
```

The result of res(6/10) in Maple 2018.1 is -0.555168575558567e-2+1.00000000000001*10^(-22)*I
with the default settings of abserr and relerr.

## Bug?...

I tried the following and lost the kernel connection.

```restart;
u:=(-5*ln(x)^4*Pi^4-20*ln(x)^2*Pi^4-8*Pi^4+120*MeijerG([[0, 0], [1, 1, 1]], [[0, 0, 0, 0, 0], []], x^Pi))/(120*Pi^4*(-1+x)^2);
evalf(u);
int(%, x = 4/10 .. 6/10, numeric);
```

I shall report that as an SCR.

## Chop...

I once had that problem and then made a small procedure I called Chop lacking a better name.
Unlike evalf[n] which will turn exact expressions like Pi or sqrt(2) into floats with n digits Chop only works on existing floats.

```Chop:=proc (x::anything, n::posint) local op1, op2, k;
if not hastype(x, float) then return x end if;
if type(x, float) then
op1, op2 := op(x);
k := length(op1);
evalf(round(op1*10^(-k+n))*10^(op2-n+k), n)
elif type(x, complex(float)) then
procname(Re(x), n)+I*procname(Im(x), n)
else
evalindets(x,float,procname,n)
end if
end proc;

```

Examples:

```Chop(123.456789,4);
Chop(1.23456789,4);
evalf[10](f^Pi);
evalf[4](%); #Actually doesn't do anything!
Chop(%,4);
```

 First 7 8 9 10 11 12 13 Last Page 9 of 150
﻿