Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



Re: Strange behavior of System Modeler
Posted:
Mar 3, 2013 10:58 PM


Hi, > > I have just started studying WSM and the Modelica language by > implementing the simple pendulum model (class) es. described on p.33 of > Peter Fritzson's book Introduction to "Modeling and Simulation of > Technical and Physical Systems with Modelica". The code (pendulum > equation written as a DAE) is: > model DAEExample "DAEExample" > constant Real PI=3.14159265358979; > parameter Real m=1,g=9.81,l=0.5; > output Real F; > output Real x(start=0.5),y(start=0); > output Real vx,vy; > equation > m*der(vx)=x/l*F; > m*der(vy)=y/l*F  m*g; > der(x)=vx; > der(y)=vy; > x^2 + y^2=l^2; > end DAEExample; > > I run the example in WSM and I get totally meaningless results. The > solver is the default one (DASSL) which, I think, is the right one > for handling DAEs. > > Any ideas?
If your meaningless results look similar to the ones I get it might be a deficiency of the WSM solver. The way it is written the problem is quite "unfriendly" for a numeric DAE solver and needs some more involved tricks to be solved, see e.g.:
<https://www.modelica.org/events/workshop2000/proceedings/old/Mattsson.pdf>
where the tricks that Dymola uses to get this solved are described. It came as a pleasant surprise to me that Mathematicas NDSolve (with the IndexReduction option) solves this correctly, so it might be a problem that WRI actually knows how to solve and you probably want to report this.
A workaround is of course to use a formulation that is more "friendly" to the solver, e.g. by reformulating in terms of polar coordinates (but I think Peter Fritzson might well use it in this form to demonstrate something...).
hth,
albert
(* Mathematica code for the above: *)
m=1;g=9.81;l=0.5; res=NDSolve[{ m*vx'[t]==x[t]/l*f[t], m*vy'[t]==y[t]/l*f[t]m*g, x'[t]==vx[t], y'[t]==vy[t], x[t]^2+y[t]^2==l^2, x[0]==0.5, vx[0]==0, y[0]==0, vy[0]==0 },{x,vx,y,vy,f},{t,0,5}, Method>{"IndexReduction">Automatic} ]
xsol=x/.res[[1]] ysol=y/.res[[1]]
Plot[{xsol[t],ysol[t]},{t,0,5}]



