Preben Alsholm

MaplePrimes Activity


These are answers submitted by Preben Alsholm

I will answer your last two questions in my preferred way:
 

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:

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.

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.

You may try depends.

depends(expr,a);

where expr is the polynomial or equation.

There could very likely be a singularity, not only a numerical problem.
I made the following changes:
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);


 

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.

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);
## Answer 7/12
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);
## Answer 1/2

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

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};

 

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

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.

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

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.

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.
 

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.

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