Preben Alsholm

13471 Reputation

22 Badges

20 years, 251 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

Eq2 is not defined in that worksheet.

@rcorless Actually odetest also works:
 

odetest(y(t)=theta,diff(y(t),t,t)+omega^2*sin(y(t))=0);  # 0

where theta is defined as in your worksheet.

@rcorless I managed to get zero by treating the two terms in the ode differently:
 

restart;
k := sin(alpha/2);
theta := 2*arcsin(k*JacobiSN(omega*t, k));
t1:=expand(omega^2*sin(theta));
t2:=simplify(diff(theta, t, t));
simplify(t1+t2);  # 0

 

@Muhammad Usman XJTU Since you say that you are new to Maple and have little knowledge of it, you probably shouldn't have just taken over the code written by somebody else for another ode.
Your ode as presented in Maple code is first order, so as pointed out by Tom Leslie and Rouben Rostamian you are only allowed one condition and shooting doesn't make any sense.
Rouben suggests that you present your problem in words and not in Maple code.
I think you should do that. Is your intended ode actually the one we see?

You don't seem to be defining the function PHIb anywhere. What is it?

Perhaps needless to say, but the user needs kernelopts(floatPi=true) as the code is written. That setting is the default.
I have chosen not to use it though,so I just replace 2.*Pi by twoPi, where twoPi is
local twoPi:=evalf(2*Pi).

@vv My guess is that he mistakenly thinks that the expression depends on the ratio L/z only.
That is not correct, which I guess is your reason for wondering.

A simple counterexample is given here:

restart;
ex := L/(z*sqrt(z^2 + L^2));
eval(ex,{L=1,z=1}); # L/z = 1
simplify(eval(ex,{L=2,z=2})); # L/z =1
## Those are not equal
##
simplify(eval(ex,L=LDZ*z)) assuming z>0; # LDZ=L/z


 

@Thomas Richard I quite agree! But in a moment of I don't know I did this:
 

restart;
eq2:=efe+x1*sinh(y1*l/r1)+x2*cosh(y1*l/r1)=0;
eq3:=efe+x1*sinh(-y1*l/r1)+x2*cosh(-y1*l/r1)=0;
solve({eq2,eq3},{x1,x2});
eq2+eq3;
eq2-eq3;

and we see that the last two lines confirm the solve result which is

{x1 = 0, x2 = -efe/cosh(y1*l/r1)}

Thus we really need to see that worksheet!

@Rouben Rostamian  As for the difficulty of using dsolve backwards in time because of instability we can use your x__ss integral instead.
For the present concrete case x__ss is even as a function of t:
 

x__ss := Int(exp(-s + 1/100*cos(10*t) - 1/100*cos(-10*t + 10*s)), s = 0 .. 10); # Using 10 instead of 7

This integral obviously can be computed numerically with equal ease for negative as well as positive t.

Picking x0 = 1.00997547504228, Digits = 15, abserr = 1e-18, relerr = 5e-13, and using 10 instead of 7 in x__ss makes the plots of x__ss and x_exact look the same on -13..10, but dsolve starts "falling down" after that. Your x__exact doesn't.

Graphing the difference between x__exact and x__ss on just t = -10 .. 2 we get


It gets much worse for on -13..2.

### Added.
Changing T to 20 we get this difference:

So using X__ss is the way to go!

@Rouben Rostamian  That is very nice indeed.
A comment:

1. It could be mentoned that while x(t) = xtr+xss is a solution to de neither xtr or xss is a solution to the de.
2. I do believe,however, that there must be a solution defined on the whole real axis. That solution is periodic. All other solutions will tend to either -infinity or +infinity as t -> -infinity.

## The argument for part of point 2 is that the the horizontal strip { (t,x) |  10/11 < x < 10/9 } is positively invariant. In fact the right hand side of the ode f(t,x) is negative for any x > 10/9 and positive for x < 10/11.
Outside the strip the solutions will grow exponentially backwards in time. The ones above the strip will go to +infinity, those below to - infinity.
Consider the x-interval  I = (10/11, 10/9). Starting at any arbitrary time t0 with x(t0) in that interval there are 3 possibilities when going backwards in time: Either the solution leaves the strip upwards to infinity (case 1), downwards to minus infinity (case 2), or stays in the strip for all times (case 3).

Using existence and uniqueness of the initial value problem for this ode we see that the solutions belonging to cases 1 and 2 correspond to non-overlapping open subintervals of the interval I.
The interval of x-values corresponding to case 3 must be closed and not empty, possibly consisting of just one point.
This still leaves open the possibilty that this last interval consisting of solutions that stay in the strip for all negative as well as positive times consists of more than one point.

In our present example surely it seems that there is only one point.


Notice that that x-point will depend on the chosen value of t0. When t0 = 0 is chosen it is something like x0 = 1.00997547504228. Because of the obvious instability in the negative direction solving numerically from (0, x0) will not give us the solution in that direction:

 

How do you define an equibrium point for a non-autonomous ode?

If you mean a constant solution, i.e. x(t) = k, where k is a constant, then it is trivial to see that your ode doesn't have such a solution.
By inserting x(t) = k into x'(t)=-(1+0.1*sin(10*t))*x(t)+1  you get 0 = -(1+0.1*sin(10*t))*k+1 = -k - 0.1*k*sin(0.1*t) +1

Since this has to hold on an interval t = a..b k must be zero. Thus we get 0 = 1, a contradiction!

What you may be able to show is that all(?) solutions tend towards some common periodic solution x0(t) as t -> infinity, i.e.
x(t)-x0(t) -> 0 as t -> infinity.

It is easy to give convincing numerical evidence of that claim.

@AHSAN I would recommend Maple's own manuals. There is a User Manual and a Programming Guide.
Go to https://www.maplesoft.com/documentation_center/

They have manuals from previous versions too.

Furthermore, I would strongly recommend that you use Maple Input (1D math input) instead of 2D and also Worksheet Mode.

To make that change in Maple go to the menu item Tools. Then Options/Display/Input display. Here choose Maple Notation.

After that and still with the Options window open go to Interface/Default format for new worksheets. Here choose Worksheet.

Now click on Apply Globally. Don't worry, the changes can easily be undone.
The changes only affects new worksheets. So try opening a new worksheet and start writing some Maple code on an input line marked with >.
###########################################

I have added a link to a worksheet in the format recommended here.

The worksheet includes the code shown above, but also conversion of the second order ode to a first order system in the usual way. Furthermore this system is solve exactly as well as numerically by a shooting method.

MaplePrimes21-08-05_ode_bvp.mw

@AHSAN You ignored my correction of your syntax error. Please see above.
Since you refer to a paper, where they use numerical solution, and not a textbook teaching you about bvps, there seems to be absolutely no point in not using the exact method given above.
I was beginning to think that your instructor had given you this problem as an exercise and wanted it done in the way you describe.

That not being so just use the result sol1 or sol (sol = sol1) as found also by Maple 2019.2.

That result is:

u(y) = -1/2*y/beta + 1/12*(4*beta*p*y + 1)^(3/2)/(beta^2*p) + 1 + 1/2*h/beta - 1/12*(4*beta*h*p + 1)^(3/2)/(beta^2*p);

From that you can compute all the values you want.
Incidentally, why didn't you update your Maple 2019.0 to 2019.2?

@vs140580 I don't have Excel on my computer, so I cannot help you with that.
I'm assuming that you know of the ExcelTools package. It has Export and Import to and from Excel.

First 14 15 16 17 18 19 20 Last Page 16 of 225