Preben Alsholm

13728 Reputation

22 Badges

20 years, 250 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@tomleslie I never used that icon, but just tried entering sqrt(1-u) after no success with sqrt(1+u).

sqrt(1-u)

No success with sqrt(1+sqrt(5)), but change the sign and voila:

sqrt(1-sqrt(5))

No success with sqrt(3+sqrt(5)), but sqrt(3-sqrt(5)) gives

sqrt(3-sqrt(5))

I shall not bother you with more examples although it is amusing to try and see what works.

What is it that is being repeated? I see absolutely nothing being repeated.

@iman Actually, my comments were intended to dissuade you from pursuing the symbolic approach and take up the numerical approach instead.

@iman Before you get started on solving the problem numerically, be sure to replace omega^2 by omega2:
sys:=subs(omega^2=omega2,{p1,p2,p3});

In addition to the suggestion from Markiyan notice that odeplot allows functions of the dependent variables too, like ln(y(t)).

res1:=dsolve({op(eval(DissMod1,Pars1))} union {y(0) = 1, B[1](0) = 1, B[2](0) = 0}, numeric);
plots:-odeplot(res1,[t,y(t)], 0..1000);
plots:-odeplot(res1,[t,ln(y(t))], 0..1000);

@Carl Love That is indeed a good question.
I tried these:
restart;
ode1:=diff(x(t),t)=1/[x(t)+1]*sin(t);
res:=dsolve({ode1,x(0)=0},numeric); #No apparent problem
plots:-odeplot(res,[t,x(t)],0..1); #Now the problem shows up
## Also of course it would show up just by doing
res(0.1);
Error, (in res) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

ode2:=diff(x(t),t)=[x(t)+1]*sin(t);
dsolve({ode2,x(0)=0},numeric); #Problem right away
## It is not surprising that these two complain right away:
dsolve({ode1,x(0)=0},numeric,range=0..1,output=piecewise);
dsolve({ode1,x(0)=0},numeric,output=Array([0.5,1]));
##
My guess would be that dsolve/numeric (with output=procedurelist/listprocedure/operator) sets up a procedure (or several) to calculate the solution, but does the minimal amount of calculation to do so. Thus it is left to chance whether an expression like [x(t)+1] (which is not syntactically wrong) is caught.

Addendum.

Setting infolevel[dsolve]:=5 shows that the computation procedure for the ode1 problem has YP[1] = sin(t)/[Y[1]+1]  (coming from DEtools[convertsys]).
That problems come up right away for ode2 is due to the result of eval(rhs(ode2),t=0), which is [0].
For ode1 we get eval(rhs(ode1),t=0), which is just plain 0.
A third version ode3:=diff(x(t),t)={x(t)+1}*sin(t);
This one also survives dsolve because eval(rhs(ode3),t=0) is also just 0.


@iman I shall elaborate on the simple example I gave above:

restart;
ode:=diff(y(x),x)=1/y(x);
dsolve(ode); #General solution
dsolve({ode,y(0)=0}); #Two solutions: i.e. lack of uniqueness.
## Trying numerically:
res:=dsolve({ode,y(0)=0},numeric);
plots:-odeplot(res,[x,y(x)],0..1); #NO LUCK (no surprise either)
#Now trying a midpoint method:
res2:=dsolve({ode,y(0)=0},numeric,range=0..1,method=bvp[midrich]); #NO LUCK
##Trying to find the parameter a so that the solution with y(1)=a satisfies y(0)=0:
res3:=dsolve({ode,y(1)=a},numeric,parameters=[a],output=listprocedure);
Y:=subs(res3,y(x));
res3(parameters=[sqrt(2)]); #According to results above a=sqrt(2) should work: CHEATING!
plots:-odeplot(res3,[x,y(x)],0..1); #Indeed it looks good.
##Now try to find a without cheating, i.e. without knowing an exact solution:
Q:=proc(a) res3(parameters=[a]); Y(0) end proc;
Q(1.4); #Error
Q(1.42); #OK
##We modify Q such that it returns -10 when an error in Y(0) occurs:
Q:=proc(a) res3(parameters=[a]);
  try
    Y(0)
  catch:  
    -10
  end try
end proc;
##This version of Q works with fsolve and fsolve finds the correct parameter:
fsolve(Q,[1..2]);
#############
The lesson to learn is that although a singularity is present in your ode with y(0)=0, there may still be a solution (here two), but it may require some extra work to find.


@iman Here is the code to find the denominators starting with your system of odes (called sys here).

restart;
sys:={diff(F(eta), eta, eta, eta)+E(eta)*(3*F(eta)*(diff(F(eta), eta, eta))-(diff(F(eta), eta))^2)/(.18*K(eta)^2)+(2*(diff(K(eta), eta))/K(eta)-(diff(E(eta), eta))/E(eta))*(diff(F(eta), eta, eta))+E(eta)/(0.9e-1*K(eta)^2) = 0,
diff(theta(eta), eta, eta)+15*E(eta)*(F(eta)*(diff(theta(eta), eta))-theta(eta)*(diff(F(eta), eta)))/K(eta)^2+(2*(diff(K(eta), eta))/K(eta)-(diff(E(eta), eta))/E(eta))*(diff(theta(eta), eta)) = 0,
diff(K(eta), eta, eta)+E(eta)*(1.5*F(eta)*(diff(K(eta), eta))-K(eta)*(diff(F(eta), eta)))/(0.9e-1*K(eta)^2)+(2*(diff(K(eta), eta))/K(eta)-(diff(E(eta), eta))/E(eta))*(diff(K(eta), eta))+(diff(F(eta), eta, eta))^2-E(eta)^2/(0.9e-1*K(eta)^2) = 0,
diff(E(eta), eta, eta)+7.2222*E(eta)*(3*F(eta)*(diff(E(eta), eta))-E(eta)*(diff(F(eta), eta)))/K(eta)^2+(2*(diff(K(eta), eta))/K(eta)-(diff(E(eta), eta))/E(eta))*(diff(E(eta), eta))+1.872*E(eta)*(diff(F(eta), eta, eta))^2/K(eta)-27.7333*E(eta)^3/K(eta)^3 = 0};
##Your boundary conditions:
bcs:={F(0) = 0, (D(F))(0) = 0, K(0) = 0, E(0) = 0, theta(0) = 1,((D@@2)(F))(1) = 0, K(1) = 0, E(1) = 0, theta(1) = 0};
indets(sys,specfunc(diff));
## The highest derivatives for F, E, K, and theta:
highD:={diff(F(eta),eta$3),diff(E(eta),eta,eta),diff(K(eta),eta,eta),diff(theta(eta),eta,eta)};
sysS:=solve(sys,highD); #Solving for the highest derivatives
##Finding the denominators (although they are in fact visible in sysS):
denom~(subs(sysS,highD));
##Answer:  {K(eta)^2*E(eta), K(eta)^3*E(eta)}
##These denominators are sure to give you numerical problems because of your boundary conditions for E and K.
## As I said, try a midpoint method and you find you still have a problem. So in fact the reported singularity at eta=0 and/or at eta=1 might indeed actually be a singularity.
#################
##Simple example:
restart;
ode:=diff(y(x),x,x)=1/y(x);
dsolve({ode,y(0)=0,y(1)=0},numeric,method=bvp[midrich]);




I'm not sure I get your point unless the equations were supposed to be differential equations, so that the first set would be:
Diff(x[1],t) = f[1](p[1]....p[n],x[1]...x[m])
...
etc.

And the second set
Diff(y[1],s) = f[1](q[1]....q[p],y[1]...y[m])
...

etc.
where also t is changed to s.
?

@vv Indeed!
The help page for sum says:

To add a finite sequence of values, rather than compute a formula, use the add command.  For example, add(k, k=0..9) returns 45.  Although the sum command can often be used to compute explicit sums, it is strongly recommended that the add command be used in programs if an explicit sum is needed, in particular, when summing over all elements of a list, Array, Matrix, or similar data structure.
(emphasis added)

@H-R Maple simplies u/sqrt(u) to u^(1/2), because it is always correct:

Proof:  Write u = r*exp(I*t), where r > 0 and  -Pi < t <= Pi.
Then  u^(1/2) = r^(1/2) * exp( I*t/2), since -Pi < -Pi/2 < t/2 <= Pi/2 <= Pi, so that t/2 is the principal argument of u^(1/2).
Thus by using the properties of the exponential function we get:

u/sqrt(u) = (r/r^(1/2))*exp(I*t -I*t/2) = r^(1/2)*exp(I*t/2) = u^(1/2).

rsolve doesn't use f(0)=1:
restart;
eq:=f(n)=f(n-1)/2+f(n+1)/2;
pb:={eq, f(0)=1, f(6)=0};
res:=rsolve(pb, {f(n)});
eval(%,n=0); #Using that f(0)=1 we get:
res2:=eval(res,f(5)=1/6);


It is not explicitly stated in the help page for rsolve, but seems to be assumed that boundary conditions are given sequentially with none missing between lowest and highest n. Since only two conditions can be given here, f(0) and f(1) could be given, or f(5) and f(6) for instance.

Example:
pb2:={eq, f(5)=1/6, f(6)=0};
rsolve(pb2,f(n));


@vv Actually, the result is not wrong, but isn't satisfactory since f(0) = 1 is not used:

eq:=f(n)=f(n-1)/2+f(n+1)/2;
pb:={eq, f(0)=1, f(6)=0};
res:=rsolve(pb, {f(n)});
eval(%,n=0); #So using that f(0)=1 we get:
res2:=eval(res,f(5)=1/6);


@vv and Markiyan: Thanks for the references.

Thank you Carl (and acer) for sharing with us this insight.

Off topic, but inspired by the fact that the false eigenvalues RPf (as remarked by Carl) seems to lie on a circle, I recalled my observation a few years ago that a randomly chosen polynomial of high degree has roots which tend to cluster around a circle of radius 1.
I wondered why, and still do.

restart;
p:=RandomTools:-Generate(polynom(integer(range=-99..99),z,degree=100));
r:=fsolve(p=0,z,complex):
plots:-complexplot([r],style=point,symbol=circle,scaling=constrained);

First 91 92 93 94 95 96 97 Last Page 93 of 230