Thank you.
Could I trouble you to give an example of your method?
this should work; it does in Maple 10 classic
>restart;
> with(plots): #I think you have to have this
> la:=1:#or whatever value you wish to have
> a:=0:# or whatever starting value you wish for fi
> b:=eval(Pi/2):#or whatever end value you wish for fi
> plot(la*cos(FI),FI=a..b);# Actually Maple will not let you use 'fi'
#the above implies that y = la*cos(FI)
In case you are not familiar with the filing syntax I have used the following
D is my drive name
MapleVariables is a folder I created in the windows browser on drive D
(I am not sure of the distinction between mpl and m files)
then I would have the folowing line
>save variable 1, variable 2, ...variable n, "D:\\MapleVariables\\somefile.m":
(I always seem to forget the last comma before "D)
and then
>read "D:\\MapleVariables\\somefile.m":
reads variable 1, variable 2, ...variable n,
I have not adjusted the solver settings apart from the absolute error which is e-4 (seems good with this) and the time step is t=dur/51 (dur<6)
The 51 is fairly arbitrary.
Using ‘vec’ to hold the time points, as above, I had set the solve procedure up with the output option,
>dsol:=dsolve(numeric,number=12,procedure=solproc,start=0,
initial=ic,procvars=dvars,method=rkf45,
abserr=Float(1,-4),output=array(vec)):
and then to store the returned for all the intermediate step values I had
> ent:=(entries(dsol))
However, my only means – within my knowledge of Maple - of getting the optimised values out the procedure was to drop the output option and call the procedure multiple times with something of the form
dsol2:=dsolve(numeric,number=12,procedure=solproc,start=0,
initial=ic,procvars=dvars,method=rkf45,abserr=Float(1,-4)):
Pmatrix:=Matrix(51,4):#
for j from 1 to 51 do
dsol2(vec[j]);
for c from 1 to 4 do
Pmatrix[j,c]:=(Pvector[c]):
end do:
end do:
Is it possible to assign the optimized values assigned to ‘Pvector’ within the procedure directly to the ‘ent’ array for
I was thinking of do a first pass to find the natural time steps and then calling those.
Thank you once again for your interest in this.
To answer my own question, perhaps an initial pass with steppast=true in order to determine the natural steps and then use the natural steps as the calculated time points?
There wouldn't be a routine in Maple that automates this?
Thank you again Allan
On my actual calculations while steppast=false allows me to get the global variables out of the procedure at the relevant time points there is a 12-fold overhead in calculation time.
Are there any obvious means of improvement? I have set abserr as high as I can.
To answer my original question, if you assign to a global variable within the numeric solving procedure using rkf45 are the values returned the values that occur for the specified final value of the independent variable? The answer would seem to be that unless you set steppast=false or happen to choose the natural step value for the independent variable the answer is, no.
I have adjusted my original file to demonstrate this.
Download 2129_for maple2.mwsView file details
Perhaps someone with a greater insight to the numerical solvers might like to confirm this.
This sheet shows a simple single body forward dynamic model in plane motion with optimisation of two applied forces to provide a smooth angular velocity and an linear acceleration constraint in one direction. There is coloumb friction force in one direction.
The model in the uploaded sheet
Download 2129_for maple.mwsView file details
isn't for any purpose. It is just to exhibit the phenomena I am trying to understand in a more cumbersome routine.
In brief I specify the duration of the simulation by the independent variable u. And of course if
dur:=2
and I have
dsol2:=dsolve(numeric,number=8,procedure=solproc,start=0,initial=ic,procvars=dvars,method=rkf45,abserr=Float(1,-4)):
and I call dsol2(dur);
the value of u on the return is u=2.
But if I save u along with the values of the forces in a global vector the value that is returned for u is not 2. It appears to always be a little more than u.
How can I then be sure that the values for the forces saved in the global vector are the ones at time u=2?
Is there something more complex about dsolve than I reckoned with.
Thank you Acer.
That now seems entirely obvious but leads to another question.
I have assigned the value of the independent variable to this extra array
along with the results of NLPSolve
But if I use eg dsolve(.9090909091e-1) the value of the independent variable that is returned is .109918045794435674. I expected the independent value that is returned to be the value with which the procedure is called.
Are there subtleties about returning the value of the independent variable or is the above an indication that the values returned by NLPSolve are for a different value of the independent variables than that used to calculated the dependent variables?
Thank you. I have learnt something!
I am using Maple 10 classic on a PC
My reference to 'machine specific' was to the fact that intermediate calculations are saved to pre-existing folders referenced in absolute terms so the files would need to be modified to run.
Thank you for your interest in this
CASE 1-Does not solve (My orginal) when run with the sub-routine as below
RTFR:=proc(vel,RT_0)
> local vn,N,k,mu,lambda,friction,v0 ;
> global RT_FR_max;
> vn:=vel:
> N:=1:
> k:=.005:
> mu:=RT_FR_max+1*10^(-20):
> v0:=0:
> if abs(vn) then v0:=(-1)*signum(vn)
> end if;
>
> friction:=-mu*N*(vn-k*v0)/(abs(vn)+k);
>
> if vn=0
> then friction:=RT_0
> #else friction:=-mu*N*(vn-k*v0)/(abs(vn)+k);
> end if;
>
> end proc:
The solver completes with
"Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up"
#############################################################################
CASE 2-Does not solve [I've added the lprint(friction)] command when run with the sub-routine as below
RTFR:=proc(vel,RT_0)
> local vn,N,k,mu,lambda,friction,v0 ;
> global RT_FR_max;
> vn:=vel:
> N:=1:
> k:=.005:
> mu:=RT_FR_max+1*10^(-20):
> v0:=0:
> if abs(vn) then v0:=(-1)*signum(vn)
> end if;
>
> friction:=-mu*N*(vn-k*v0)/(abs(vn)+k);
>
>
> if vn=0
> then friction:=RT_0
> #else friction:=-mu*N*(vn-k*v0)/(abs(vn)+k);
> end if;
> lprint(friction);
>
> end proc:
The solver terminates with
"
-.25
-.25
-.25
-.25
Error, (in f) unable to store '.99882335063022e-1+(NULL)' when datatype=float[8]
Error, invalid input: entries received dsol, which is not valid for its 1st argument, t"
#############################################################################
CASE 3-Does solve when run with the sub-routine as below
RTFR:=proc(vel,RT_0)
> local vn,N,k,mu,lambda,friction,v0 ;
> global RT_FR_max;
> vn:=vel:
> N:=1:
> k:=.005:
> mu:=RT_FR_max+1*10^(-20):
> v0:=0:
> if abs(vn) then v0:=(-1)*signum(vn)
> end if;
>
> friction:=-mu*N*(vn-k*v0)/(abs(vn)+k);
>
>
> if vn=0
> then friction:=RT_0
> #else friction:=-mu*N*(vn-k*v0)/(abs(vn)+k);
> end if;
> lprint(friction);
> friction:=friction;
> end proc:
.25000000000000000001
-.22754962511063422927
-.22752734255366036904
-.13013463640892057384
-.13022861982265447511
-.18198071247644153013
-.18194930791954201510
-.25000000000000000001
The solver behaves correctly and the results are coherant with the physical problem
The numerical values are hundreds of lines long the above are just the last lines. The correct
magnitude for RT_0 is 0.25 and the sign is correct.
You really need the files don't you?
There are a number of files working together which are set up to write to specific locations on my machine so I won't trouble you with the actual files.
Here is a sample of the code. As you will see below I have modified a subprocedure called from the main integration procedure (solproc) and this works fine. I would be interested to understand why.
#THIS WORKS#
RLLS:=proc(vel,RL_0)
local vn,N,k,mu,lambda,friction,v0 ;
global RL_LS_max;
vn:=vel:
N:=1:
k:=.005:
mu:=RL_LS_max+1*10^(-20):
v0:=0:
if abs(vn)
I am using
dsol:=dsolve(numeric,number=12,procedure=solproc,start=0,initial=ic,procvars=dvars,method=rkf45,abserr=Float(1,-4),output=array(vec)):
It is a solution for equations of motion.
With my limited knowledging of error finding I'm using
printlevel:=1000;
to examine why the routine does not solve.
There in the thousands of lines of output from printlevel I find repeatedly
ERROR "unable to store %1 when datatype=%2"
I assumed that somehow I was returning a variable with the
wrong datatype to the main routine.
More generally is there any information on how to interpret what printlevel outputs?
Thanks for the interest
1 2 3 4 |
Page 3 of 4 |