Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Kamel Algeria In your example the result is the same if you don't use assumptions.

In the help page for solve it says (emphasis added):

"If the output of the solve command is a piecewise-defined expression, then the assuming command can be used to isolate the desired solution(s). If the output is not piecewise-defined, in particular, if the output is constant, assumptions on the independent variables may be ignored. If there are parameters in the input equations, the solve command will use those assumptions in its computations."

It would help tremendously if you would provide us with the lines of code you have used, so we don't stumble around in the dark.

@goli As I mentioned earlier, you could turn your equation into a differential equation. That I have done at the end of this comment.

If you use fsolve, then add a starting value to ensure that you get the positive branch. I have chosen H=1 (it shouldn't be too crucial, which positive value you choose, but don't make it too small).

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
Y := z->fsolve(eq(z), H=1):
plot(Y, -2.2 .. 5);
implicitplot(eq(z), z = -3 .. 10,H=-10..10,gridrefine=2);
plot(fdiff(Y, [1], z), z = -2.2 .. 6, caption = "The derivative of y using fdiff on Y");
yp := implicitdiff(eq(z), H, z);
plot(eval(yp, H = 'Y(z)'), z = -2.2 .. 6, caption = "The derivative of y using implicitdiff and Y");
plot(eval((1+z)*yp/H, H = 'Y(z)'), z = -1 .. 5);

#The differential equation approach:

ode:=diff(H(z),z)=subs(H=H(z),yp);
eval(eq(0),H=1);
res:=dsolve({ode,H(0)=1},numeric,output=listprocedure);
odeplot(res,[z,H(z)],-2.22..5);
Ynum:=subs(res,H(z));
Ynum(2.345678);


@goli As I mentioned earlier, you could turn your equation into a differential equation. That I have done at the end of this comment.

If you use fsolve, then add a starting value to ensure that you get the positive branch. I have chosen H=1 (it shouldn't be too crucial, which positive value you choose, but don't make it too small).

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
Y := z->fsolve(eq(z), H=1):
plot(Y, -2.2 .. 5);
implicitplot(eq(z), z = -3 .. 10,H=-10..10,gridrefine=2);
plot(fdiff(Y, [1], z), z = -2.2 .. 6, caption = "The derivative of y using fdiff on Y");
yp := implicitdiff(eq(z), H, z);
plot(eval(yp, H = 'Y(z)'), z = -2.2 .. 6, caption = "The derivative of y using implicitdiff and Y");
plot(eval((1+z)*yp/H, H = 'Y(z)'), z = -1 .. 5);

#The differential equation approach:

ode:=diff(H(z),z)=subs(H=H(z),yp);
eval(eq(0),H=1);
res:=dsolve({ode,H(0)=1},numeric,output=listprocedure);
odeplot(res,[z,H(z)],-2.22..5);
Ynum:=subs(res,H(z));
Ynum(2.345678);


@goli Instead of me guessing what you are doing, I suggest that you bring the whole code beginning with a restart to indicate that nothing has been defined before.

@goli Instead of me guessing what you are doing, I suggest that you bring the whole code beginning with a restart to indicate that nothing has been defined before.

I cannot reproduce your result of 858.000000000000114.

I get the integer 858 (not 858.) as I expected. I'm using Maple 14, worksheet interface, Maple input (of course).

If I change the division by the integer 1000 to division by the float 1000. I get 858.000000000000.

@goli Have you tried it yourself? Did you run into any problems? If so, what were they?

@goli Have you tried it yourself? Did you run into any problems? If so, what were they?

@goli Surely the positive solution is used also for y'(x).

Try this plot, which shows the function and its derivative.

plot(['Y(x)',eval(yp,y='Y(x)')],x=-3..1.5,caption="y (i.e. Y, red) and the derivative of y (using implicitdiff and Y, blue)",color=[red,blue]);

Notice that

factor(yp);

returns 5*y*(4*x+7)*(1+x)^2/(9*y^2+2+7*x+9*x^2+5*x^3+x^4).

Thus yp has two zeros, nicely confirmed by the plot.

Another way of getting to your numerical solution is to turn your equation into a differential equation with an appropriate initial condition, which could be y(-1) = 1, since x=-1 and y=1 satisfies eq(x).

ode:=diff(y(x),x)=subs(y=y(x),yp);

res:=dsolve({ode,y(-1)=1},numeric,output=listprocedure);
Yode:=subs(res,y(x));
plot(Yode,-3..1.5);

If you only want to plot y and y', then you need only do

plots:-odeplot(res,[[x,y(x)],[x,rhs(ode)]],-3..1.5);

@goli Surely the positive solution is used also for y'(x).

Try this plot, which shows the function and its derivative.

plot(['Y(x)',eval(yp,y='Y(x)')],x=-3..1.5,caption="y (i.e. Y, red) and the derivative of y (using implicitdiff and Y, blue)",color=[red,blue]);

Notice that

factor(yp);

returns 5*y*(4*x+7)*(1+x)^2/(9*y^2+2+7*x+9*x^2+5*x^3+x^4).

Thus yp has two zeros, nicely confirmed by the plot.

Another way of getting to your numerical solution is to turn your equation into a differential equation with an appropriate initial condition, which could be y(-1) = 1, since x=-1 and y=1 satisfies eq(x).

ode:=diff(y(x),x)=subs(y=y(x),yp);

res:=dsolve({ode,y(-1)=1},numeric,output=listprocedure);
Yode:=subs(res,y(x));
plot(Yode,-3..1.5);

If you only want to plot y and y', then you need only do

plots:-odeplot(res,[[x,y(x)],[x,rhs(ode)]],-3..1.5);

@goli I think you are using Robert Israels's reformulation of the equation

eq2:= x->(y^2 - (1+x)^3 - (1+x)^4)^5 - y=0;

which for y>0 is equivalent to eq. However, for general y they are not equivalent.

The following does give two curves, yes, one above the x-axis, and one below.

plots:-implicitplot(eq2(x), x=-10..10,y=-10..10,grid=[100,100],gridrefine=2);

Using eq you only get a curve above the x-axis.

I thought you were interested in the one above?

@goli I think you are using Robert Israels's reformulation of the equation

eq2:= x->(y^2 - (1+x)^3 - (1+x)^4)^5 - y=0;

which for y>0 is equivalent to eq. However, for general y they are not equivalent.

The following does give two curves, yes, one above the x-axis, and one below.

plots:-implicitplot(eq2(x), x=-10..10,y=-10..10,grid=[100,100],gridrefine=2);

Using eq you only get a curve above the x-axis.

I thought you were interested in the one above?

@goli Here are the different commands.

restart;
eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;
Y:=x->fsolve(eq(x),y);
plot(Y,-10..10);
plots:-implicitplot(eq(x), x=-10..10,y=0..120,grid=[150,150]);
plot(fdiff(Y,[1],x),x=-2..2,caption="The derivative of y using fdiff on Y");
yp:=implicitdiff(eq(x),y,x);
plot(eval(yp,y='Y(x)'),x=-2..2,caption="The derivative of y using implicitdiff and Y");

@goli Here are the different commands.

restart;
eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;
Y:=x->fsolve(eq(x),y);
plot(Y,-10..10);
plots:-implicitplot(eq(x), x=-10..10,y=0..120,grid=[150,150]);
plot(fdiff(Y,[1],x),x=-2..2,caption="The derivative of y using fdiff on Y");
yp:=implicitdiff(eq(x),y,x);
plot(eval(yp,y='Y(x)'),x=-2..2,caption="The derivative of y using implicitdiff and Y");

First 216 217 218 219 220 221 222 Last Page 218 of 229