Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

Did you try computing the integral before attempting minimization?

int(0.1e-3+.5*t+0.2e-2*t^2-b*t-a, t = 0 .. 300);

Output:

@acer I tried changing initial condition in your version as I did in my comment above. It doesn't work.

Just checking existing conditions doesn't work either:

ONE_SOLUTION(initial);

The initial condition for D(h) can be included in the parameters, however:

ODE_SOLUTION := dsolve({ODE, h(0) = 0,D(h)(0)=y1}, numeric, method=rosenbrock_dae,parameters = [y1,A, B, C, E, F]):
ODE_SOLUTION('parameters'=[Pi/2,0,0,1,1,1]):






@ There is ambiguity in determining the initial condition even in the case you consider.
The initial condition for u, i.e. u(0)=u1, satisfies cos(u1)=0. Thus there is a choice to be made.
After having found the solution as given, we can do

sol(initial=[0,0,-Pi/2]);

and when plotting we find that we have another solution.


@Carl Love I couldn't agree more!

@Alejandro Jakubi Thank you very much, I shall try to contact Gaston Gonnet.

@Markiyan Hirnyk Yes, but you are using two calls to eval.
You cannot do
eval(f(a*x), x = a, a = 9);

@acer Using method=nonlinearsimplex worked in Maple 15 on my simple example, which is of the same kind as the one described by joha.

@acer I noticed that joha uses Maple 14. On a simple example I tested, this device doesn't work in Maple 15 (and thus most likely not in Maple 14). It does work in Maple 18 (and in Maple 16).

I think it might be necessary to use the operator form with a user supplied gradient as shown by you in 2010 (July 14).

@Benrath plot3d(x*y,x=y..1,y=0..1);

@nm I think you are being rather unreasonable.
When I do
?pi
it brings up a help page, which starts out saying
"For information on 3.14159..., see Pi."

I think that that is as helpful as it gets!

And that
?foo
brings up a help page about Units,length, (foot)
just shows that the help system is trying to be helpful.

In your worksheet you don't attempt any solution of the system (sys), but you are finding equilibria by using solve. I don't get any error at all, but the output:



i.e. infinitely many equilibria, here parametrized by v, x, and z.

If you want to solve sys, you will undoubtably have to do it numerically. For that the constants m and q have to be made concrete and initial (or boundary) values must be known.

@sharena2 Why do you think results will get better using shooting?

However, if you insist, I suggest shooting from eta=L, which will only require two parameters.
For comparison I continue from the code I gave above. Thus 'sol' is the output from dsolve/bvp.

#Shooting from the right with the 2 parameters fL and vL:
ICSL:={F(L) = 0, H(L) = 2, f(L) = fL, g(L) =-fL, u(L) = 0,v(L)=vL};
resL:=dsolve(ODE union ICSL,numeric,parameters=[fL,vL],output=listprocedure);
#We shall need these two:
ff,uu:=op(subs(resL,[f,u](eta)));
#Getting parameter values from sol (no shooting):
pp:=subs(sol(L),[f(eta),v(eta)]);
#Quick graphical test:
resL(parameters=pp):
plots:-display(Array([[seq(plots:-odeplot(resL,A[i],0..L),i=1..3)],[seq(plots:-odeplot(resL,A[i],0..L),i=4..6)]]));

#There are 2 parameters  (fL and vL) and 2 requirements: f(0) = 0 and u(0) = 1
pL:=proc(fL,vL) option remember;
  resL(parameters=[fL,vL]);
  [ff(0),uu(0)-1]
end proc:
p1:=proc(fL,vL) pL(_passed)[1] end proc:
p2:=proc(fL,vL) pL(_passed)[2] end proc:
#Testing pL on the parameters from sol (pp):
pL(op(pp)); #Should be close to zero
ppL:=fsolve([p1,p2],pp);
pp;
ppL-pp;
#With no so good a guess:
fsolve([p1,p2],[.4,0]);
resL(parameters=%):
plots:-display(Array([[seq(plots:-odeplot(resL,A[i],0..L),i=1..3)],[seq(plots:-odeplot(sol,A[i],0..L),i=4..6)]]));


You normally would do something like (example)

if not type(aaa,numeric) then error "%1 should be of type numeric",aaa end if;

Do you have any reason whatsoever to think that finding an analytical solution is possible?

Using phi and phi0 instead of your 2D versions:

eq1:=diff(phi(tau), tau)-1-a(tau)^2*cos(2*phi(tau))-a(tau)^2+cos(phi0-phi(tau))/a(tau)-cos(2*phi(tau)) = 0;
eq2:=diff(a(tau), tau)-a(tau)^3*sin(2*phi(tau))+a(tau)-sin(phi0-phi(tau))-a(tau)*sin(2*phi(tau)) = 0;
#dsolve({eq1,eq2}); #Doesn't seem to end
#Just try setting a=1 to see the complexity of this:
eq3:=eval(eq1,a=1);
dsolve(eq3);
value(%);

@mapleus The problem lies with fsolve, or rather with the difficulty of finding an accurate numerical solution to something like x^3 = 0. The derivative at the solution is zero. Thus Newton's method, which relies on the derivative at the approximate point, will get into trouble. Surely, fsolve is set up to handle that kind of thing, but apparently fails to find a solution that it itself is satisfied with.
Here is an example obtained by taking n:=3 and i:=1.

restart;
n:=3: i:=1:
u:=(x1,x2,z)->(x1+x2*z+9*z^i)^n;
p:=proc(x1,x2) option remember; local res;
  res:=[u(x1,x2,1),u(x1,x2,2)];
  #print([_passed,res]); It will print an awful lot if you remove #
  res
end proc:
p1:=proc(x1,x2) p(_passed)[1] end proc;
p2:=proc(x1,x2) p(_passed)[2] end proc;
fsolve([p1,p2],[0,0]); #Responds quickly unevaluated
T:=op(4,eval(p)): #The remember table of p (huge)
#eval(T); #To see it remove #
L:=entries(T,nolist):
AL:=map(x->abs(x[1])+abs(x[2]),[L]): #The sum of absolute values of the error in making zero
min(AL); #The minimal error
##Or you could sort AL
SAL:=sort(AL):
SAL[1..5];
nops(AL);
You see that fsolve is rather reluctant to declare success!
It is not too difficult to find the input corresponding to the minimal error (i.e. the best point (x1,x2) ), but since this is only an illustration I stop here.

First 133 134 135 136 137 138 139 Last Page 135 of 231