Preben Alsholm

13471 Reputation

22 Badges

20 years, 298 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@necron When I download from the link I still get t1.
However, I changed t1 to t myself and had no luck: In fact the system lost contact with the kernel when I tried your worksheet.

You are trying to evaluate the double integral numerically, but what is t1?

@red1eco You need concrete values for the parameters. I chose them randomly as positive floating point numbers between 1 and 9. You also need initial values. I took them out of my head:
restart;
sys2 := [diff(x(t), t) = r*x(t)*(K-x(t)-a*y(t))/K, diff(y(t), t) = s*y(t)*(L-y(t)+b*x(t))/L-q*y(t)*E(t), diff(E(t), t) = W*((p-tau)*l*y(t)-c)*E(t)];
params := `minus`(indets(sys2, name), {t});
r:=rand(1..9.0);
r(); #Test of the random number generator r
paramsR:={seq(p=r(),p=params)};
#Your sys2 is a list. By sys2[] you get the sequence of elements in the list:
res:=dsolve(eval({sys2[],x(0)=1.2,y(0)=0.5,E(0)=.1},paramsR),numeric);
plots:-odeplot(res,[t,x(t)],0..1); #Just the graph of x
plots:-odeplot(res,[[t,x(t)],[t,y(t)],[t,E(t)]],0..1,thickness=3,legend=[x,y,E]); #All in the same plot


@Markiyan Hirnyk Markiyan, I'm not about to enter that old discussion with you about stability of equilibria for nonlinear autonomous systems. These results are extremely well known by people who have had a decent course in differential equations. Since I know that you are Ukrainian, I believe that you speak Russian (I don't). There are quite a few very good textbooks on odes written by Russians. I know that because they have also been translated into English.

Why not use LinearAlgebra:-Eigenvector?

Choosing y ' (0) = 0 as suggested by tomleslie works fine numerically on the interval 0..1:

restart;
edo := -b*y(x)*(diff(y(x), x))^2*x-a*y(x)*x+(diff(y(x), x, x))*x+2*(diff(y(x), x)) = 0
edo1 := subs(a = 1, b = 1, edo);
res:=dsolve({edo1,y(0)=1,D(y)(0)=0},numeric,method=bvp[midrich],range=0..1);
plots:-odeplot(res,[x,y(x)],0..1);

## To push it further, you have to work harder:

res:=dsolve({edo1,y(0)=1,D(y)(0)=0},numeric,method=bvp[midrich],range=0..2,approxsoln=[y(x)=1-0.3*x^2],initmesh=256);
plots:-odeplot(res,[x,y(x)],0..2);

The solution to the singular initial value problem {edo1, y(0)=1, D(y)(0)=0} probably tends to infinity at some finite x0, i.e. y(x)-> infinity for x-> x0 from the left. x0 seems to be greater than 2 and probably less than 3.

Note. There is some evidence for an analytic solution at zero having radius of convergence less than 2.141580437. Powers all even. Using the terms up to degree 50 of that solution (sol) as approxsoln in dsolve you can get success with
res:=dsolve({edo1,y(0)=1,D(y)(0)=0},numeric,method=bvp[midrich],range=0..2.09,approxsoln=[sol],initmesh=1024);
plots:-odeplot(res);






@ecterrab We agree that a syntax like f(a)(b) certainly is useful as in e.g. D(f)(b).
I was deliberately vague when calling f(t)(t) a weird construction.

But notice that it was the appearance of several "weird" constructions like that that led me to find Torleif's mistake, which he admits to having made, when making last minute changes to a worksheet.

So far nobody asked you why you need this function. What do you intend to use it for? The answer to that may determine which method is preferable.

@ecterrab Sure, but can you give a non-weird example of f(t)(t)? Point being that there are 2 t's.

D(sin)(sin) would not make much sense, would it?

@torabi What sense does it have to remove nonlinear terms?

@torabi I had a nagging suspicion that your change to b=0 was not a good idea: If the system has the solution s=0, g=0, no what the values of omega11 and omega22 are,  then using b=0 would most likely lead to that solution.

I didn't get down to checking. But now by doing

plots:-odeplot(res[indx[i]],[[x,s(x)],[x,g(x)]],0..1,thickness=3);

for i = 1, 2, 3, ..., 11 (one at a time), you will see that 1,6,7,9,10, and maybe 11, are probably (s,g)=(0,0) with numerical roundoff.

So I suggest changing b=~0 to b=~1.

@torabi You already have them in your latest worksheet where you used the code



If you want a tag put on them, you can do this instead::

seq(subs(res[i](0),[i,omega11(0),omega22(0)]),i=indx);

Your problem is nonlinear, so I don't see any possibility for an analytical solution to the problem.

@torleian You say "I'm sorry, but I don't quite understand your earlier objection".
All my objections have been about the point that something like f(t)(t) is strange, although not syntactically wrong in Maple.

@torleian Let me state first that I'm not a user of the Physics package.
I find the objectionable term in g0. I still find it meaningless, to put it bluntly.

So to me it is no wonder that the system must object at some point. Your procedure diff_dI has to differentiate g0 w.r.t. theta(t) (and also w.r.t. diff(theta(t),t)) using the Physics version of diff.

I tried the two loops below. You get the error messages we have seen, but are the results (where there is no error) as you expected?
I find all of these expressions in IDTS below strange:
(cos(y(t)+phi(theta(t))))(t), (diff(theta(t), t))(t), (diff(theta(t), t, t))(t), (diff(y(t), t))(t), (diff(y(t), t, t))(t), (sin(y(t)+phi(theta(t))))(t), ((D@@2)(phi))(theta(t)), ((D(phi))(theta(t)))(t), (((D@@2)(phi))(theta(t)))(t).

Simple example: What does (diff(y(t), t))(t) mean?

g0;
IDTS:=indets(g0,function);
for ttt in IDTS do
  try
    res:=diff(ttt,theta(t))
  catch:
    print(lasterror);
  end try;
  Diff(ttt,theta(t))=res
end do;

for ttt in IDTS do
  try
    res:=diff(ttt,diff(theta(t),t))
  catch:
    print(lasterror);
  end try;
  Diff(ttt,diff(theta(t),t))=res
end do;

###################
This reminds me of an error some of my students would make occasionally, but which didn't seem to have consequences when plotting:

restart;
f:=sin(t); #f is not defined as a procedure (function) so should not be used as such.
plot(f,t=0..Pi); #Correct syntax
plot(f(t),t=0..Pi); # Incorrect, but works!
## So why does it work? Because numbers are accepted as constant functions:
5(7); #output 5
f(6);
subs(t=5,f(t));
%;




It seems that you want to find the jacobian matrix for the right hand sides of the system ode1,ode2 below.
I would do it like this:

restart;
f1:=(x,y)->x^2+y^2;
f2:=(x,y)->x^2-y^2;
ode1:=diff(x(t),t)=f1(x(t),y(t));
ode2:=diff(y(t),t)=f2(x(t),y(t));
VectorCalculus:-Jacobian([f1(x,y),f2(x,y)],[x,y]);

First 92 93 94 95 96 97 98 Last Page 94 of 225