Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@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.
 

@nm I just downloaded the Physics update from the Cloud and in Maple 2018.2 64bit on a computer with Windows 10 and last updated a day ago.
I experienced no problem during download. The pde problems posted by Edgardo Cheb-Terrab recently on MaplePrimes were solved with no errors.

@ecterrab Yes I didn't mention it, but did indeed try it. What I hoped for was that Maple picked the correct branch in each case itself.
Actually, I only needed splitting at x=Pi/2. But as it turned out that would have hidden the use of the wrong branch on Pi/2..Pi since the two integrals would be zero both of them regardless of branch.

@Rouben Rostamian  You are right. Finding the inverse of cos restricted to an interval different from 0..Pi seems not to be implemented in Change, which uses solve.
solve gives the same result for all four:
 

solve({cos(2*x)=z,x>0,x<Pi/4},x);
solve({cos(2*x)=z,x>Pi/4,x<Pi/2},x);
solve({cos(2*x)=z,x>Pi/2,x<3*Pi/4},x);
solve({cos(2*x)=z,x>3*Pi/4,x<Pi},x);

It has been a long time since I used Mathematica, but I remember it issuing a warning about having to invert transcendental functions.

@Mariusz Iwaniuk I'm not convinced. The method used is ftoc (try infolevel[int]:=5), but the antiderivative is discontinuous.
No indication that that is taken into account.

u:=convert(cos(2*x)/(1+2*sin(3*x)^2), exp);
U:=int(u,x);
plot(Re(U),x=-1..Pi+1);


 

@vv If we use the integral in the form int(-cos(2*x)/(-2+cos(6*x)), x ) instead, it appears by plotting (!) that that antiderivative is continuous. As pointed out in my answer Maple doesn't find the definite integral with the integrand written like that though.

@Carl Love Thanks, Carl, for your detailed explanation. I think I got it now. I'm slow.

@ecterrab The other day I deleted the Physics Update folder (named 2018) and after that the warning disappeared.
The point of that exercise was to find out if the warning would go away. Since it did, the warning will only come up for people actually downloading the Physics update. Thus they would presumably know that they have such a thing.
Thus I'm not too worried anymore.

Clearly, the hack above can be turned into a procedure.
Input the dsolve solution (res) and the number of odes of the corresponding first order system (ndiff).

restart;
FreezeParams:=proc(res,ndiff::posint,$) local R,params,ind,i,j,S;
  params:=rhs~(res('parameters'));
  ind:=[ListTools:-SearchAll('undefined',params)];
  R:=convert(eval(res),string);
  S:=NULL;
  j:=0;
  for i from 1 to nops(params) do
    if not member(i,ind) then 
      j:=j+1;
      S:=S,cat("Y",[ndiff+j])=convert(params[i],string)
    end if;
  end do;
  R:=StringTools:-Subs({S},R);
  parse(R)
end proc:
###  
## Example.
res:= dsolve({diff(y(x),x,x)=a*y(x)+b, y(0)=c,D(y)(0)=d}, numeric, parameters= [a,b,c,d]):
res0:=FreezeParams(res,2); 
res(parameters);
res0(parameters);
res(parameters=[a=1]);
res1:=FreezeParams(res,2); 
res0(parameters);
res1(parameters);
res(parameters=[b=2,c=3]); # Two more set
res2:=FreezeParams(res,2);
res2(parameters);
res2(parameters=[7,8,9,10]); #Attempting to set parameters
res2(parameters); # Not set
res0(parameters); #Unchanged
res1(parameters); # Unchanged
res(parameters=[d=4]); # All set.
res3:=FreezeParams(res,2);
result:=res(1); 
res3(1);
res(parameters=[0,0,0,0]);
res(1);
res3(1); # same as 'result'.

 

@Carl Love It seems that the solution procedure in dsolve/numeric uses Y[2] for the parameter in the case of one parameter in a system of one equation:
 

restart;

res:= dsolve({diff(y(x),x)=a*y(x), y(0)=1}, numeric, parameters= [a]):
res(parameters=[2]);
result:=res(1);
R:=convert(eval(res),string):
R2:=StringTools:-Subs("Y[2]"="2.",R):
res2:=parse(R2):
res2(1); # As result above
res2(parameters);  # Query
res2(parameters=[7]); # Seemingly set
res2(1); # But it isn't

Guessing the (obvious) generalisation:

restart;
### First order, two parameters:
res:= dsolve({diff(y(x),x)=a*y(x)+b, y(0)=1}, numeric, parameters= [a,b]):
res(parameters=[1,2]);
result:=res(1);
R:=convert(eval(res),string):
R2:=StringTools:-Subs({"Y[2]"="1.","Y[3]"="2."},R):
res2:=parse(R2):
res2(1);
res2(parameters); 
res2(parameters=[7,9]);
res2(1);
######################################
restart;
### Second order, two parameters:
res:= dsolve({diff(y(x),x,x)=a*y(x)+b, y(0)=1,D(y)(0)=0}, numeric, parameters= [a,b]):
res(parameters=[1,2]);
result:=res(1);
R:=convert(eval(res),string):
R2:=StringTools:-Subs({"Y[3]"="1.","Y[4]"="2."},R):
res2:=parse(R2):
res2(1);
res2(parameters); 
res2(parameters=[7,9]);
res2(1);


 

First 37 38 39 40 41 42 43 Last Page 39 of 231