<

Preben Alsholm

12813 Reputation

22 Badges

17 years, 310 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@mary120 In this version I have shown how to find acceptable initial values for given values of delta.
The clue is that ode has diff(n(x),x)^2 which is nonnegative. Thus the rest of lhs(ode) must be negative (or zero) since the right hand side of ode is 0.

restart;  #NEW V

Digits:=15:
V:=-n^2*((T[e]+T[i])/T[e])*((((n)*(1+1/(2*delta[p]))))-ln(n)-ln(n/(2*delta[p]))-1-1/(4*delta[p])-delta[p]);
W:=eval(V,{T[e]=6/5,T[i]=1/100,n=n(x)}); # Omit delta[p]
ode:=diff(n(x),x)^2+2*W=0;
############################

############## When does a real solution exist for ode?
E:=lhs(eval(ode,{n(x)=n})); # E must be negative or zero for a real solution of ode
sol:=solve(E,n,allsolutions);
E1:=eval(sol,_Z2=0); # If you run the first time after restart you see _Z2. If you repeat without restart you get _Z4
E2:=eval(sol,_Z2=-1);
mMa:=eval([E1,E2],delta[p]=0.7);
mMb:=eval([E1,E2],delta[p]=0.9)
plot([E1,E2],delta[p]=0.7..0.9);
plot(eval(E,delta[p]=.7),n=0..2);
plot(eval(E,delta[p]=.9),n=0..2);
mMa;
mMb;
m:=min(mMa[1],mMb[1]);
M:=max(mMa[2],mMb[2]);
##So the initial values acceptable for ode and ODE if delta[p] values are between 0.7 and 0.9
##must be chosen outside of the values in mMa and nMb since E must be negative (or zero). 
##So n(0) must lie in the open interval 0..m or in the open interval M..infinity.
############################
# Example of an acceptable start valid for delta[p] between 0.7 and 0.9:
eval(convert(ode,D),{n(x)=m-.01}); # Start at 0 < n < m
ICS:=solve(eval(%,x=0),{D(n)(0)});
ic:=n(0)=m-.01;
ODE:=diff(ode,x);
## Using the parameters option to dsolve/numeric:
RESpar:=dsolve({ODE,op(ICS[1]),ic},numeric,parameters=[delta[p]]);

RESpar(parameters=[delta[p]=0.7]); # # Setting the parameter
plots:-odeplot(RESpar,[x,n(x)],-5..5); p1:=%: 
RESpar(parameters=[0.9]);# You can omit the name of the parameter when setting it.  
plots:-odeplot(RESpar,[x,n(x)],-5..5,color=blue); p2:=%:
plots:-display(p1,p2);
## A procedure that helps doing the above:
Q:=proc(dp,{range::range:=-10..10}) if not dp::realcons then return 'procname(_passed)' end if;
   RESpar(parameters=[dp]); # Setting the parameter
   plots:-odeplot(RESpar,[x,n(x)],range,_rest)
end proc;
# Simple input:
Q(0.7);
Q(0.9,color=black,range=-5..5,linestyle=3);
## An animation in delta[p]:
interface(warnlevel=0); # To avoid all those warnings.
plots:-animate(Q,[delta[p],range=-5..5],delta[p]=0.7..0.9);
plots:-animate(Q,[delta[p],range=-5..5,color="DarkGreen"],delta[p]=0.7..0.9,trace=24);

@mary120 You can do as I have posted earlier. If you insist on n(0) = 1 you cannot consider delta[p]>=1/2:
 

restart;  #NEW V

Digits:=15:
V:=-n^2*((T[e]+T[i])/T[e])*((((n)*(1+1/(2*delta[p]))))-ln(n)-ln(n/(2*delta[p]))-1-1/(4*delta[p])-delta[p]);
W:=eval(V,{T[e]=6/5,T[i]=1/100,n=n(x)}); # Omit delta[p]
ode:=diff(n(x),x)^2+2*W=0;
############################

eval(W,delta[p]=0.7);
############################
eval(convert(ode,D),{n(x)=1});
ICS:=solve(eval(%,x=0),{D(n)(0)});
ic:=n(0)=1;
ODE:=diff(ode,x);
eval(ODE,delta[p]=0.7);
## Using the parameters option to dsolve/numeric:
RESpar:=dsolve({ODE,op(ICS[1]),ic},numeric,parameters=[delta[p]]);

RESpar(parameters=[delta[p]=1/5]); # # Setting the parameter
plots:-odeplot(RESpar,[x,n(x)],-3..2,view=0..20); p1:=%: 
RESpar(parameters=[1/4]);# You can omit the name of the parameter when setting it.  
plots:-odeplot(RESpar,[x,n(x)],-3..2,view=0..20,color=blue); p2:=%:
plots:-display(p1,p2);
## A procedure that helps doing the above:
Q:=proc(dp,{range::range:=-10..10}) if not dp::realcons then return 'procname(_passed)' end if;
   RESpar(parameters=[dp]); # Setting the parameter
   plots:-odeplot(RESpar,[x,n(x)],range,_rest)
end proc;
# Simple input:
Q(1/5);
Q(1/4,color=black,range=-3..2,view=0..20,linestyle=3);
## An animation in delta[p]:
interface(warnlevel=0);
plots:-animate(Q,[delta[p],range=-5..5,view=0..20],delta[p]=0.1..0.49);
plots:-animate(Q,[delta[p],range=-5..5,view=0..20,color="DarkGreen"],delta[p]=0.1..0.49,trace=24);

What is Delta?  Even if you think it's obvious, tell us anyway, please.

@mary120 You cannot use n(infinity) = 0 ever. You must use a finite value for x0 (and n0) in n(x0) = n0.
In your case you cannot even have n(x0) = 0 because ln(n(x)) appears in your ode.

The introduction of the ode for x, which was indeed odd, turns out to be superfluous (a euphemism for stupid).
So just get rid of that and replace X(x) with x.
Here is the revised code:

restart;

Digits:=15:
V:=n^2*((T[e]+T[i])/T[e])*((((n-1)*(1+1/(2*delta[p]))))-ln(n)-ln(n/(2*delta[p])));
W:=eval(V,{T[e]=6/5,T[i]=1/100,n=n(x)}); # Omit delta[p]
ode:=diff(n(x),x)^2+2*W=0;
############################
eval(convert(ode,D),{n(x)=1});
ICS:=solve(eval(%,x=0),{D(n)(0)});
ic:=n(0)=1;
ODE:=diff(ode,x);
## Using the parameters option to dsolve/numeric:
RESpar:=dsolve({ODE,op(ICS[1]),ic},numeric,parameters=[delta[p]],output=listprocedure);
###
N:=subs(RESpar,n(x));
RESpar(parameters=[delta[p]=1/5]); # # Setting the parameter
plots:-odeplot(RESpar,[x,n(x)],-1..2);
nm:=numapprox:-infnorm(N,0..2,'x00');
plots:-odeplot(RESpar,[x-x00,n(x)],-1..5);
####
Qtrans:=proc(dp,{range::range:=-5..5,scene::list:=[x-x0,n(x)],normrange::range:=-2..2,print::truefalse:=true}) local x00;
if not dp::realcons then return 'procname(_passed)' end if;
   RESpar(parameters=[dp]); # Setting the parameter
   numapprox:-infnorm(N,normrange,x00); #
   if print then printf("Translated by %f\n ",x00) end if; 
   plots:-odeplot(RESpar,eval(scene,x0=x00),range,_rest)
end proc;
###########
Qtrans(0.49);

Qtrans(0.49);

Qtrans(0.49,normrange=-0.1..0.1);# The normrange is not large enough, so x00 is too small 

Qtrans(1/10);

Qtrans(1/10,normrange=0..2);

plots:-animate(Qtrans,[delta[p]],delta[p]=0.1..0.49);

plots:-animate(Qtrans,[delta[p],color="DarkGreen",print=false],delta[p]=0.1..0.49,trace=24);

The results are the same as in "An attempt" above.

@mary120 You can achieve a translated version of the last code by modifying as follows.
The basic idea is to introduce an ode for x (! yes indeed).
This time the animated trace version looks like this:

restart;

Digits:=15:
V:=n^2*((T[e]+T[i])/T[e])*((((n-1)*(1+1/(2*delta[p]))))-ln(n)-ln(n/(2*delta[p])));
W:=eval(V,{T[e]=6/5,T[i]=1/100,n=n(x)}); # Omit delta[p]
ode:=diff(n(x),x)^2+2*W=0;
############################
eval(convert(ode,D),{n(x)=1});
ICS:=solve(eval(%,x=0),{D(n)(0)});
ic:=n(0)=1;
ODE:=diff(ode,x);
odeX:=diff(X(x),x)=1;
icX:=X(0)=0;
## Using the parameters option to dsolve/numeric:
RESpar:=dsolve({ODE,odeX,op(ICS[1]),ic,icX},numeric,parameters=[delta[p]],output=listprocedure);
###
N,XX:=op(subs(RESpar,[n(x),X(x)]));
RESpar(parameters=[delta[p]=1/5]); # # Setting the parameter
plots:-odeplot(RESpar,[x,n(x)],-1..2);
nm:=numapprox:-infnorm(N,0..2,'x00');
plots:-odeplot(RESpar,[X(x)-x00,n(x)],-1..5);
Qtrans:=proc(dp,{range::range:=-5..5,scene::list:=[X(x)-x0,n(x)],normrange::range:=-2..2,print::truefalse:=true}) local x00;
if not dp::realcons then return 'procname(_passed)' end if;
   RESpar(parameters=[dp]); # Setting the parameter
   numapprox:-infnorm(N,normrange,x00); #
   if print then printf("Translated by %f\n ",x00) end if; 
   plots:-odeplot(RESpar,eval(scene,x0=x00),range,_rest)
end proc;
Qtrans(0.49);
Qtrans(0.49,normrange=-0.1..0.1);# The normrange is not large enough, so x00 is too small 
Qtrans(1/10);
Qtrans(1/10,normrange=0..2);
plots:-animate(Qtrans,[delta[p]],delta[p]=0.1..0.49);
plots:-animate(Qtrans,[delta[p],color="DarkGreen",print=false],delta[p]=0.1..0.49,trace=24);

@mary120
Comments on tha question:
All solutions will be periodic, but not with the same period.
So to get the full picture for each curve we only need to plot over one period.
Your initial condition n(0) = 1 will have to be given up if you want all the maximal to occur at x=0.
#########
Furthermore, remember that the maxima occur at points where the ode doesn't satisfy requirements for uniqueness of solutions. So attempts at starting at a point where a maximum occurs can easily run into numerical problems (and seems to do so).
You could say that my second approach to your ode just happens to work the way you want it.
But there are in fact infinitely many solutions passing through a given maximum:
Take a solution on the way to the top, reaching the maximum, deciding to stay there for a while and then go down to the minimum and then go up etc.

@mary120 Here is how you can do that:
 

restart;

Digits:=15:
V:=n^2*((T[e]+T[i])/T[e])*((((n-1)*(1+1/(2*delta[p]))))-ln(n)-ln(n/(2*delta[p])));
W:=eval(V,{T[e]=6/5,T[i]=1/100,n=n(x)}); # Omit delta[p]
ode:=diff(n(x),x)^2+2*W=0;
############################

############################
eval(convert(ode,D),{n(x)=1});
ICS:=solve(eval(%,x=0),{D(n)(0)});
ic:=n(0)=1;
ODE:=diff(ode,x);
## Using the parameters option to dsolve/numeric:
RESpar:=dsolve({ODE,op(ICS[1]),ic},numeric,parameters=[delta[p]]);

RESpar(parameters=[delta[p]=1/5]); # # Setting the parameter
plots:-odeplot(RESpar,[x,n(x)],-10..10); p1:=%: 
RESpar(parameters=[1/4]);# You can omit the name of the parameter when setting it.  
plots:-odeplot(RESpar,[x,n(x)],-10..10,color=blue); p2:=%:
plots:-display(p1,p2);
## A procedure that helps doing the above:
Q:=proc(dp,{range::range:=-10..10}) if not dp::realcons then return 'procname(_passed)' end if;
   RESpar(parameters=[dp]); # Setting the parameter
   plots:-odeplot(RESpar,[x,n(x)],range,_rest)
end proc;
# Simple input:
Q(1/5);
Q(1/4,color=black,range=-5..5,linestyle=3);
## An animation in delta[p]:
plots:-animate(Q,[delta[p],range=-5..5],delta[p]=0.1..0.49);
plots:-animate(Q,[delta[p],range=-5..5,color="DarkGreen"],delta[p]=0.1..0.49,trace=24);

Here is the first plot displaying p1 and p2:

And here is the animation with trace:

@Joe Riel I suppose that what TechnicalSupport wanted to define was:
 

find_vars_in_proc:=proc(f :: procedure, $)
  return {op(2, eval(f))};
end proc;

 

@Anthrazit That line from the help page 'Initially Known Names' is the exact same as in Maple 12 and probably goes back even longer. OK a comma in Maple 12 has been changed to a period later.
The statement is:
" Pi     math constant pi (pi); use Pi for calculations.   
           evalf(Pi) is approximately 3.14159265... "

I think the point is that the lower case pi is used first, and then it has to be pointed out that Pi is the actual name of the mathematical constant in Maple.

@tomleslie I don't think that conversion to radians in Maple 2022 takes place in general (and it didn't in Maple 2021 either).
Maple 2022 appears to have a problem when using Units:-Simple:-`<` where one of the arguments is zero.
The code for Units:-Simple:-`<` was different in Maple 2021 and worked fine.

 

with(Units[Simple]);
alpha:=45*Unit(degree);
beta:=30*Unit(degree);
evalb(alpha>beta); # true
evalb(Pi/4>Pi/6);  #Not decided as expected
evalb(alpha>0); #  0<Pi/4 in Maple 2022, but true in Maple 2021

I would consider this a bug in Units:-Simple:-`<` in Maple 2022.

@Anthrazit Apparently min and max work, but not `<`.

restart;
with(Units[Simple]);
alpha:=45*Unit(degree);
if alpha=0 then "alpha=0" elif max(0,alpha) = alpha then "alpha>0" end if;

So you could define LessThan like this:

LessThan:=proc(x,y)
       if x-y=0 then return false end if;
       if min(x-y,0)=x-y then true else false end if
end proc;

@harry4939 Save your worksheet. Then use the icon with fat green arrow pointing up in the MaplePrimes editor  to upload the worksheet.
I uploaded my worksheet Test.mw, which was made in Maple 2022.0. No problems.
Test.mw

 

An approach that doesn't begin by solving for the derivative, but differentiates the ode and thereby turns it into a second order ode:
 

restart;

Digits:=15:
V:=n^2*((T[e]+T[i])/T[e])*((((n-1)*(1+1/(2*delta[p]))))-ln(n)-ln(n/(2*delta[p])));
W:=eval(V,{T[e]=6/5,T[i]=1/100,delta[p]=1/5,n=n(x)});
ode:=diff(n(x),x)^2+2*W=0;
############################
eval(convert(ode,D),{n(x)=1});
ICS:=solve(eval(%,x=0),{D(n)(0)});
ic:=n(0)=1;
ODE:=diff(ode,x);
RES1:=dsolve({ODE,op(ICS[1]),ic},numeric);
plots:-odeplot(RES1,[x,n(x)],-10..10);
RES2:=dsolve({ODE,op(ICS[2]),n(0)=1},numeric);
plots:-odeplot(RES2,[x,n(x)],-10..10);
#############################

The first plot:

You forgot to attach the worksheet showing the result from Maple 11.

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