Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Zeineb Could you upload a worksheet in which you have tried to implement my code?
Then I can have a look at it.

@Zeineb As Tom Leslie points out the code I gave works in Maple 18.02 as well. Elementwise operations were introduced in Maple 13. So the code should work in all versions from 13 and up.
Since you report the error message
" Error, (in dsolve/numeric/process_input) input system must be an ODE system, got independent variables {t, u} "

clearly you didn't just copy and paste the code into a fresh worksheet. The 'u' in "independent variables {t, u}"  must come from the integration variable used in your integrals.
The integrals should have disappeared totally in SYS.

restart;
dsolve(diff(y(x),x)=x);

Result:

  Warning, this version of the "Physics Updates" is intended

   for Maple release 2018.1, while you are using Maple 2018.2

  . Please check at http://www.maplesoft.com/products/maple/feat\

  ures/physicsresearch.aspx for the availability of a version

   of the "Physics Updates" for the copy of Maple you are using.
                              1  2      
                       y(x) = - x  + _C1
                              2         

The warning would be confusing to many. What has physics got to do with this ode one might ask.
I know the answer,  but I still find the warning unfortunate.

A consolation is that the warning only appears once after a restart. The next (or same) dsolve command doesn't result in a warning.

Followup:
I just downloaded Physics Updates from the Cloud (yes I know a link was given in the first warning).
The warning this time was:
"`Warning, this version of the "Physics Updates" is intended for Maple release ``2018.2`` and _libraryversion <= `1356656`, while you are using Maple release ``2018.2`` with _libraryversion = `1362973`. `

So maybe the earlier warning should have included a warning not to use the Cloud :-)

Note: I have now tried the link and as I thought downloaded and installed the Physics update:
Result: The same warning as from the Cloud.

@Lottie The first line in the procedure Q serves the purpose of making Q(A) return unevaluated if A is not an actual real constant (or doesn't evaluate to one).

if not A::realcons then return 'procname(_passed)' end if;

You could have written this instead:
 

if not type(A,realcons) then return 'Q(A)' end if;

If the name ttxx is unassigned then Q(ttxx) will return (literally) Q(ttxx).
The reason I wrote the line the way I did is simply that it's the way I usually do it.
It allows for inclusion of all arguments actually passed, but that is irrelevant here.

@Lottie That you could do to assure yourself of your results, but remember that your Lyapunov function V is just one that works, i.e. proves local asymptotic stability of (0,0).
The region in the plane where dV is negative doesn't necessarily correspond to the basin of attraction. In this case it certainly doesn't.

While I can't explain what I see from the lack of context, I hazard a guess that it may have to do with the fact that series remembers.
 

@ERA The OP wrote:
" I am working on an assignment and have some doubts about my solution to plotting the position of an anharmonic oscilliator with the Runge-Kutta method. "
Thus he might find it comforting that dsolve/numeric seems to agree with his results.

@ERA Surely you are right, but the OP said:
" I am working on an assignment and have some doubts about my solution to plotting the position of an anharmonic oscilliator with the Runge-Kutta method. "
Thus I thought it might be a comfort to him that dsolve/numeric appears to agree.

@Christian Wolinski I have no current access to Maple 5, but the behavior in Maple 8 is the same as it is in Maple 2018.

You really should use LinearAlgebra:-MatrixExponential(A,t) instead of DETools:-matrixDE .
The latter uses the very old matrix structure dating back to before Maple 6 (!).

@sand15 For A,B, and C try:
 

addressof(eval(A));
addressof(eval(B)); # As above for A
addressof(eval(C)); # Different

So C evaluates to a different copy of the table than A and B, which presumably is the purpose of 'copy'.

@Carl Love In the worksheet the contribution from phi disappears because the parameter a is set to zero in the loop.
The parameter a is taken from
 

Va := [0, 0];

The worksheet executes without any problem and as presented on my computer with Maple 2018.1.
## It works the same in Maple 2017.3 and also Maple 12.02.
## I'm using the 64 bit version on Windows 10.

@Carl Love A comment about results from dsolve/numeric/bvp spurred by your reply above.
By using the approximate solution
 

approxsoln=[f(eta)=lambda,diff(f(eta),eta)=eta/8-1,diff(f(eta),eta,eta)=0,g(eta)=0,theta(eta)=1]

I got results in the special case that were extremely good compared to the known exact solution in the special case K=0.
Using the same approximate solution and with abserr=1e-12 I got for Table 3, column 1:
 

[2.99473305260864, 3.07341413295920, 3.20904849825216, 3.38805933695604, 3.59894177889025]

The values reported in the paper are:

[2.98754,3.06345,3.19543,3.36981,3.57500]

Thus the differences are:
 

[0.7e-2, 0.10e-1, 0.14e-1, 0.18e-1, 0.24e-1]

So I agree with you that the numbers given in that column are not accurate to 4 decimal places.

## The approximate solution used I got from turning the ode for f into a system of first order in the special case. I noticed that that made dsolve give the desired solution without a user provided approxsoln. The initial profile generated by dsolve/numeric/bvp was in the special case:
 

[f(eta) = 3, f1(eta) = (1/8)*eta-1, f2(eta) = 0]

when %infinity=8, lambda=3 and the conversion was diff(f(eta),t)=f1(eta), diff(f1(eta),eta)=f2(eta).

@666 basha OK, so why not do that yourself?

@Carl Love Let us assume that F as presented indeed solves the ode for f with the given values of the parameters, in particular with lambda = 3. Then the given boundary condition f(0) = -lambda is WRONG and must be replaced by f(0) = lambda.
For that problem there are 2 solutions for f (at least 2). The first is found directly by dsolve/numeric, the second is found by dsolve/numeric with approxsoln=[f(eta)=F]:

restart:
Digits:= 15:
ODEs:= [
   #Eq 9/18/24:
   (1+K)*diff(f(eta),eta$3) + f(eta)*diff(f(eta),eta$2) - 2*n/(n+1)*diff(f(eta),eta)^2 -
   2/(n+1)*M*diff(f(eta),eta) + K*diff(g(eta),eta) + 2/(n+1)*sigma*theta(eta)  =  0,
   #Eq 10/19/25:
   (1+K/2)*diff(g(eta),eta$2) + f(eta)*diff(g(eta),eta) - 
   (3*n-1)/(n+1)*diff(f(eta),eta)*g(eta) - 2/(n+1)*K*(2*g(eta)+diff(f(eta),eta$2))  =  0,
   #Eq 11/20/26:
   diff(theta(eta),eta$2) + Pr*f(eta)*diff(theta(eta),eta)  =  0
]:   
params:={K= 0, Pr= 0.733, sigma= 0, lambda= 3, c= 1, M= 1, n= 1,%infinity=8};
odef:=eval(ODEs[1],params);
bcsf:=f(0) = lambda, D(f)(0) = -1, D(f)(%infinity)=0;
#Eq 32:
F:= simplify(eval[recurse](lambda - (1-exp(-eta*z))/z, [z= (lambda+sqrt(lambda^2+4*M-4))/2, lambda= 3, M= 1]));
odetest(f(eta)=F,{odef,bcsf}); # OK when lambda=3
res1:=dsolve(eval({odef, bcsf},params),numeric);
plots:-odeplot(res1,[eta,diff(f(eta),eta)]); p1:=%:
res2:=dsolve(eval({odef, bcsf},params),numeric,approxsoln=[f(eta)=F]);
plots:-odeplot(res2,[eta,diff(f(eta),eta)],color=blue); p2:=%:
plots:-odeplot(res2,[eta,diff(f(eta),eta)-diff(F,eta)]);
plots:-display(p1,p2);

The two solutions for diff(f(eta),eta):

The existence of two solutions could maybe be considered a consequence of the numerically necessary replacement of infinity by a finite number (here 8).
Notice that f(eta) = lambda - eta also satisfies odef and f(0) = lambda, D(f)(0) = -1, but not D(f)(infinity) = 0.
The acceptable solution of the two presented above must be the one for which f drops right away so that f ' (eta) appears as tending to zero as eta - > infinity.

4 5 6 7 8 9 10 Last Page 6 of 196