Preben Alsholm

13728 Reputation

22 Badges

20 years, 242 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

@hirnyk If you use 1d-input there is no problem. I just copied Jean-Jacques' lines without any problem into a fresh worksheet with 1d-input.

After replacing zz1= ..  by zz1:= ..  you can do

plot([seq(zz1[1..-1,[1,k]],k=2..4)]);

or

plot([seq(zz1[1..-1,[1,k]],k=2..4)],style=point, symbolsize=20,symbol=[circle,diamond,cross]);




Your system

ex1:={diff(x(t), t, t)-1/2*x(t)*diff(x(t), t)^2 = -5, diff(theta(t), t, t)+2*diff(theta(t), t)*diff(x(t), t)/x(t) = 0};

has a problem at t=0 since x(0)=0 appears in the denominator in the ode for theta.

You can try with a very small value for x(0):

ic := {theta(0)=0, x(0)=0.001,  D(theta)(0)=0, D(x)(0)=0};
dsol := dsolve(ex1 union ic, numeric, range=-1..1);
dsol(.5);
plots:-odeplot(dsol,[t,x(t)],-1..1);

However, notice that theta(t) = 0 for all t clearly is a solution together with x(t) determined from the first equation:

icx := {x(0)=0,  D(x)(0)=0};
dsolx := dsolve({ex1[1]} union icx, numeric, range=-1..1);

plots:-odeplot(dsolx,[t,x(t)],-1..1);


Notice that isolate only works on one equation or one expression (which is interpreted as expression = 0).

Examples:

eq1:=x+a*y=b;
eq2:=2*x+y=c;
isolate(eq1,y);
solve(eq1,{y});
op(%);

# Use of elementwise operation to isolate x in both equations:
isolate~({eq1,eq2},x);

# Use of elementwise operation to isolate y in eq1 and x in eq2:

isolate~([eq1,eq2],[y,x]);

#Solving the system

solve({eq1,eq2},{x,y});

You could do like this:

restart;
ww:= [seq( [n, add(1/(n+1), k=ceil(1/3*n)..floor(2/3*n))], n=1..60)]:
plot(ww, style=point):
plot(1/3,0..60,color=blue):
plots:-display(%,%%);

Your definition of q[1,i] in the loop doesn't make it a function, i.e. procedure.

So change the lines to:

restart;
q[1,0]:=x->x^2;
for i from 1 to 2 do
q[1,i]:= 2* D(q[1,i-1]);
od;

Piecewise defined functions are not yet implemented apparently. As Georgios Kokovidis shows in his answer you can consider g on each interval of continuity or use the Optimization package. 

Here is a somewhat programmed version of the former:

g:=x->(piecewise(0 <= x, 7.5, 0))-(piecewise(0.5 < x,9, 0))-(piecewise(0.4 < x and x <= 0.6, 30*x-12, 0.6 < x,6, 0));
h:=simplify(g(x));
h:=convert(h,rational);
L:=sort(convert(discont(h,x),list));
L1:=[op(L),1];

eps:=.000001:
seq(maximize(g(x),x=L1[i]+eps..L1[i+1]-eps),i=1..nops(L1)-1);
max(%),min(%);

You can start by making the LinearAlgebra package available:

with(LinearAlgebra):

so that Maple now knows Determinant (= your det).

@MapleFans001 If you really insist that Diff(f3(t), t$2) should appear instead of Diff(f4(t), t$2) as I asked above, then certainly you have a problem with the initial conditions: The system doesn't contain Diff(f4(t), t$2).

Going directly to DEtools[convertsys] you get an error message, which with your system is not too surprising:

DEtools[convertsys](ex1,ic,[f1(t),f2(t),f3(t),f4(t)],t,y,yp);
Error, (in DEtools/convertsys) unable to convert to an explicit first-order system

However, clearly the trivial solution f1(t) = f2(t) = f3(t) = f4(t) = 0 (for all t) is a solution.

In the help page for alias we find "you cannot define one alias in terms of another".

But if you want to use alias you can do as follows:

restart;
alias(T = T(x, y),k=k(T(x,y)));
diff(T,x);
diff(k, x);

You may also want to take a look at PDEtools[declare], which works differently.

My first answer about numerical solution used the default output = procedurelist. odeplot doesn't care.

However, often I find it useful to use output = listprocedure, in which case the output is a list of procedures, and not just one procedure which outputs a list.

resLP:=dsolve({eval(motion,{k=.6,g=10}),f(0)=0,D(f)(0)=0},numeric,output=listprocedure);

#To get the procedures for f and its derivative df/dtheta you can do like this:

F,Fp:=op(subs(resLP,[f(theta),diff(f(theta),theta)]));

F(5.4321);
Fp(5.4321);

#The graph of f (i.e. F):

plot(F,0..2*Pi);
# This graph is also produced by odeplot without first fishing out F:

plots:-odeplot(resLP,[theta,f(theta)], 0..2*Pi);

# If it isn't the graph of f you want you can surely plot whatever you like using either just odeplot or F.

# If you want to change the values of k and/or g you can use the parameters option:

respar:=dsolve({motion,f(0)=0,D(f)(0)=0},numeric,parameters=[k,g],output=listprocedure);

#Now set the parameters to whatever you like:
respar(parameters=[.8,9]);
# and plot:
plots:-odeplot(respar,[theta,f(theta)],0..2*Pi);

#If you like to animate, you can define a simple procedure p like this:
p:=proc(kk,gg) global respar,g,k,f,theta;
if not type([kk,gg],list(numeric)) then return 'procname(_passed)' end if;
respar(parameters=[kk,gg]);
plots:-odeplot(respar,[theta,f(theta)],0..2*Pi)
end proc:

plots:-animate(p,[k,9.8],k=.1..0.8);
plots:-animate(p,[.5,g],g=1..10);


Well, actually dsolve produces an answer, but it is very complicated, so no answer if you ask for a solution satisfying initial conditions.

You can try solving numerically.

#Exact general solution:

dsolve(motion);

#With initial conditions it is not surprising that you don't get a result:

dsolve({motion,f(0)=0,D(f)(0)=0});

#Numerical solution (with k = 1, g = 10:

resnum:=dsolve({eval(motion,{k=1,g=10}),f(0)=0,D(f)(0)=0},numeric);
plots:-odeplot(resnum,[theta,f(theta)],0..2*Pi);

Use unevaluation quotes to prevent Maple from complaining that m is not an integer (which of course it is not, it is a letter!):

sum('x(m)',m=1..3);

or better, use add, which has special evaluation rules (like seq and mul) and is made for summations where the limits are in fact integers (like here 1 and 3):

add(x(m),m=1..3);

Maple most likely does what I suggested doing manually: Considers dx/dr and then the zero appearing in the equation for dr/dx is a problem.

We replace 0 by epsilon:

FOI := diff(r(x), x) = piecewise(r(x) < 0, 1, r(x) > 0, epsilon);
dsolve({FOI, r(0) = -1});

eval(%,epsilon=0);

In this particular case you can solve the ode for dx/dr instead (sort of):

eq2:=diff(x(r),r)=piecewise(r < 0, 1,1/epsilon);
resx:=dsolve({eq2,x(-1)=0});
resr:=solve(eval(resx,x(r)=x),r);
resr0:=eval(resr,epsilon=0);
plot(resr0,x=0..2,thickness=2);

Edited: Thanks to hirnyk for catching a typo. It should be correct now.

First 141 142 143 144 145 146 147 Last Page 143 of 160