Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Carl Love Even in document mode and with 2d input I get -1 as I should.

@Jodelkoenig I just tried in Maple 2017.3 with restart separated as Carl describes. I had no problems.
Which Maple 2017 do you have?
Mine is:
interface(version);
`Standard Worksheet Interface, Maple 2017.3, Windows 10, September 27 2017 Build ID 1265877`
and

kernelopts(version);
`Maple 2017.3, X86 64 WINDOWS, Sep 27 2017, Build ID 1265877`

@awass There is no file attached.

@vv Nice, clear, and short.

I think that in the line
 

P := indets( f, name ) minus { x }:

'minus {x}' should be removed. My reason is that the line defining Q:
 

# Solve for unknowns in g.
       Q := indets( g, name ) minus P:

is intended to get the unknown parameters in g and thus shouldn't include x.
## Actually you can replace P and Q by:

Q:={seq(a[k],k=0..n/2),r};


The question: You are allowing for several solutions. Does it ever occur?
 

@Josci95 You wrote that the code was supplied by Maple support. From the use of foldl, and the keen use of unevaluation quotes in things like 'parameters' = [ ':-a' = a, ':-b' = b ] I would have guessed that myself.
All very well, but why use foldl instead of plain add?

The worst part though, is that the code appears to be terribly inefficient.
In the uploaded worksheet below I have compared their approach to my traditional approach otherwise sticking to their settings. I used 1000 time values from the time interval 0..5. In the original there were only 5 values.
The essential difference between their approach and mine is that I have only one procedure Q. That procedure uses listprocedure output from dsolve.

solN := dsolve( { de, ic }, 'numeric','abserr'=1e-15, 'maxfun'=0, 'parameters'=[a,b],output=listprocedure ):
Xp,Yp:=op(subs(solN,[x(t),y(t)]));
###
Q:=proc(a,b) local j,ssq;
  if not [a,b]::list(realcons) then return 'procname(_passed)' end if;
  try
    solN(parameters=[a,b]);
    ssq:=add( (Xp(T[j])-X[j])^2,j=1..N)+add( (Yp(T[j])-Y[j])^2,j=1..N);
  catch:
    ssq:=10^15
  end try;
  ssq
end proc:

MaplePrimes19-01-08_ode_fit_test.mw
 

@Josci95 You could ask if there is any way of knowing that x(t) depends on the product a*b only without solving exactly as done here. In general that could be difficult if not impossible. But in this case we can see it by eliminating y(t) from the system.
Since we are pretending that we don't have data for y(t) that isn't a bad idea.
 

restart;
de := diff( x(t), t ) = a * y(t), diff( y(t), t ) = b * x(t);
ic := x(0) = 1, y(0) = 0;
###
diff(de[1],t);
ode_x:=subs(de[2]*a,%); ## depends on a*b only.
ics_x:=x(0)=1,D(x)(0)=0:  #The latter follows from ic and de[1].

 

@Rariusz Since theta(t) and diff(theta(t),t) are dependent on the groups c1, c2, and c3 only, the values of J, b, K, R, and L cannot be determined uniquely from data involving only theta(t) and diff(theta(t),t).
You will need data on i(t) too.

Example of the nonuniqueness: You can set R and L to whatever you want then J,b, and K are determined. Keep in mind though that the derivation of the 3rd order ode assumed that J, K, and L were nonzero.
 

## c-results:
Cr:= [c1 = 0.07085248546343231, c2 = 1072.44657764141, c3 = 91326.7791941768];
## sbs from my code:
sbs:={J = (L^2*c2 - L*R*c1 + R^2)/(L^3*c3^2), K = (L^2*c2 - L*R*c1 + R^2)/(c3*L^2),
    b = (L^2*c2 - L*R*c1 + R^2)*(L*c1 - R)/(L^4*c3^2)};
## Now pick L<>0 and R as you please:
eval(sbs,{R=1,L=2,Cr[]});
## or
eval(sbs,{R=10,L=7,Cr[]});

If you want you can even make b=0 for any L<>0 if you pick R = L*c1:
 

eval[recurse](sbs,{R=L*c1,Cr[]});

@Rariusz You write:
" Column 1 - time, column 2 - input signal, column 3 - position of rotor, column 4 - speed of rotor ".

Just for the record:
My interpretation related to eq1, eq2 is:

"input signal" =   V(t), "position of rotor" = theta(t), "speed of rotor"  = diff(theta(t),t).

@Rariusz Before you wanted to determine J,b,K,L,and R. Now they are gone.
I'm lost!

@Rariusz I have looked at the data. I'm not sure what are in the columns.
Right now I shall assume that the first column is time = t, the second V(t), the last to i(t) and theta(t), but in which order I wouldn't care to guess.
Before looking closer I should like to know if I'm interpreting the first 2 columns correctly and also what the last 2 are.
You mention as an example that V(10) = 11.345 , but that doesn't agree with the data if t and V correspond to the two first columns, so maybe you just took that out of your head for purposes of explanation?
### The code below concerning fitting V data has been changed.

restart;
eq1:=J*diff(theta(t),t,t)+b*diff(theta(t),t)=K*i(t);
eq2:=L*diff(i(t),t)+R*i(t)=V-K*diff(theta(t),t);
ICs:= theta(0)=0, D(theta)(0) = 0, D(theta)(0) = 0, i(0) = 0;
## Getting your data from my computer:
M:=ImportMatrix("F:/MapleDiverse/servo_open_ident_1_A.csv",source=csv);
op(1,M); # reports 5001,4
## Having a look at the 10 first rows:
M[1..10,..];
## Having a look at the 10 last rows:
M[-10..-1,..];
##Plotting what I assume to be V(t):
plot(M[..,1..2],labels=[t,V]); 
## Model below has been edited:
## Fitting those data to the model V(t) = a+b*sin(c*t+d):
V1:=Statistics:-NonlinearFit(a+b*sin(c*t+d),M[..,1..2],t,initialvalues={a=0,b=0.5,c=1/50,d=0});
## The plot of V1 is pretty close to plot(M[..,1..2],labels=[t,V]); 
## Thus V1 could be used as a simple replacement for V(t)
## at least for a trial run.
plot(M[..,[1,3]]); # plot of the third column as a function of time
plot(M[..,[1,4]],size=[1800,default]);  # plot of the fourth column as a function of time

The fit is not all that bad with the exception of t=0:
 

w1:=unapply(V1,t)~(M[..,1]):
w:=M[..,2]:
dw:=w-w1:
eval(V1,t=0);
dw[1];
(max,min)(dw[2..]); # Leaving out the difference at t=0: 0.00082,-0.00023

 

 

The third column:

The fourth column:

@Rariusz You write:
" For example in csv file I have: time vector, input vector and response of real system. I want to use this data and ode equation for parameter estimation. "

So you also have on file the response of a real system besides the 2 vectors.

Now my question is:  Which parameters are to be determined from the data?

@Rariusz OK, but I'm puzzled by your statement that you have a time vector and input vector V in a file.

Does that mean that to each element in the time vector there corresponds an element in the vector V?
So that V in your ode system should be considered a function of time rather than a constant?
If so, what role does the ode system play here?

@Ronan I agree. Once a genuine question has been answered, generally it shouldn't be possible to delete it.
By answered I just mean having received a response that deals materially with the subject matter, which of course should be related to Maple.
The OP has put the question up for people to think about and at least after having received the response mentioned he doesn't own the question anymore.
 

@nm I see that for you Physics:-Version() reports the location of Physics in the Maple 2018 lib folder.
That should just mean that the updates are not installed. They ought to be in the equivalent place reported here for my updates:
"C:\Users\Bruger\maple\toolbox\2018\Physics Updates\lib\Physics Updates.maple", `2018, December 21, 13:9 hours, version in the MapleCloud: 264, version installed in this computer: 264`
where 'Bruger' is my user name.
Notice that Tom Leslie reports the equivalent location of his Physics Updates.
That the message you get reports "... version installed in this computer: 264"  is puzzling, but may just be one thing that "worked" during your attempts at installation. But why does it say: "... version in the Maple Cloud: unable to determine" ?
Surely you know all this, so sorry: not much help.
 

2 3 4 5 6 7 8 Last Page 4 of 196