> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> in message
> <1599729044.54502.1263282254826.JavaMail.root@gallium.
> mathforum.org>...
> > > I have a set of differential equations to solve.
> I
> > > have created an ode solver compatible function in
> > > order to solve these equations, this is a nested
> > > function in order to access several variables and
> is
> > > shown below:
> > >
> > > IC2 = [0 0 0];
> > >
> > > [T2 Y2] = ode15s(@feaode, Tspan2, IC2, options);
> > >
> > > function ydot = feaode(t, y)
> > > % feaode: nested function for solving the dynamic
> > > amic snapper
> > > % system
> > > %
> > > % y = [ xA;
> > > % vA;
> > > % I ]
> > >
> > > J = y(3,1) ./ conductorArea;
> > >
> > > % interpolate position and velocity from
> input
> > > xTtemp = interp1(times,xT,t);
> > > vTtemp = interp1(times,vT,t);
> > >
> > > % relative position necessary to determine
> > > ermine forces and flux linkage
> > > xR = xTtemp  y(1,1);
> > >
> > > remxRVTaup = rem((xR/Taup),1);
> > >
> > > % now check if pole number is even or odd
> > > if round(rem((xR/Taup)remxRVTaup, 2)) == 0
> > > polysign = 1;
> > > else
> > > polysign = 1;
> > > end
> > >
> > > % We will get the numerical derivative of the
> > > of the flux linkage from the
> > > % polynomial using numerical
> differentiation
> > > psidot = polysign *
> > > sign * psidotfun(abs(rem(xR,Taup)./Taup),
> abs(J),
> > > p_psi);
> > >
> > > % determine the forces due to the magnets and
> > > ts and electrical forces at
> > > % the relative position xR with the current
> > > urrent values of J
> > > Ffea = polysign * polyvaln(p_FEAFy,
> > > FEAFy, [abs(rem(xR,Taup)./Taup), abs(J)]);
> > >
> > > % determine the spring forces on the armature
> > > mature at position xA
> > > Fs = Fs_Snapper(y(1,1), ks);
> > >
> > > % determine the emf in the coils
> > > EMF = (y(2,1)  vTtemp) .* psidot;
> > >
> > > % find the derivative of the coil current
> > > Idot = (EMF  y(3,1).*R) ./ L;
> > >
> > > ydot(1,1) = y(2);
> > > ydot(2,1) = (Ffea + Fs) ./ massA;
> > > ydot(3,1) = Idot;
> > >
> > > end
> > >
> > > This works fine, but I would like to extract the
> > > values in the variable EMF in the above function
> at
> > > each time step for examination alongside the
> values
> > > in the 'y' vector. How exactly would I need to
> > > reformulate this to do that? I have floundered
> about
> > > with mass matrices (which I don't really know
> > > anything about, but saw in another post here as
> the
> > > solution) and tried the following:
> > >
> > > % Added another initial condition, now solving
> for
> > > four variables
> > > IC2 = [0 0 0 0];
> > >
> > > % Added mass matrix to options
> > > options = odeset('Mass', [1 0 0 0; 0 1 0 0; 0 0 1
> 0;
> > > 0 0 0 0], ...
> > > 'MStateDependence','none',
> > > MStateDependence','none', ...
> > > 'MassSingular', 'yes');
> > >
> > > [T2 Y2] = ode15s(@feaode, Tspan2, IC2, options);
> > >
> > > function ydot = feaode(t, y)
> > > % feaode: nested function for solving the dynamic
> > > amic snapper
> > > % system
> > > %
> > > % y = [ xA;
> > > % vA;
> > > % I;
> > > % EMF ]
> > >
> > > J = y(3,1) ./ conductorArea;
> > >
> > > % interpolate position and velocity from
> input
> > > xTtemp = interp1(times,xT,t);
> > > vTtemp = interp1(times,vT,t);
> > >
> > > % relative position necessary to determine
> > > ermine forces and flux linkage
> > > xR = xTtemp  y(1,1);
> > >
> > > remxRVTaup = rem((xR/Taup),1);
> > >
> > > % now check if pole number is even or odd
> > > if round(rem((xR/Taup)remxRVTaup, 2)) == 0
> > > polysign = 1;
> > > else
> > > polysign = 1;
> > > end
> > >
> > > % We will get the numerical derivative of the
> > > of the flux linkage from the
> > > % polynomial using numerical
> differentiation
> > > psidot = polysign *
> > > sign * psidotfun(abs(rem(xR,Taup)./Taup),
> abs(J),
> > > p_psi);
> > >
> > > % determine the forces due to the magnets and
> > > ts and electrical forces at
> > > % the relative position xR with the current
> > > urrent values of J
> > > Ffea = polysign * polyvaln(p_FEAFy,
> > > FEAFy, [abs(rem(xR,Taup)./Taup), abs(J)]);
> > >
> > > % determine the spring forces on the armature
> > > mature at position xA
> > > Fs = Fs_Snapper(y(1,1), ks);
> > >
> > > % determine the emf in the coils
> > > EMF = (y(2,1)  vTtemp) .* psidot;
> > >
> > > % find the derivative of the coil current
> > > Idot = (EMF  y(3,1).*R) ./ L;
> > >
> > > ydot(1,1) = y(2);
> > > ydot(2,1) = (Ffea + Fs) ./ massA;
> > > ydot(3,1) = Idot;
> > >
> > > % *********** Additional Code
> *************
> > > ydot(4,1) = y(4,1)  EMF
> > >
> > > end
> > >
> > > I did this based on what I saw here:
> > >
> > >
> http://www.mathworks.nl/matlabcentral/newsreader/view_
> > > thread/246590
> > >
> > > Running the new code however returns the
> following
> > > error:
> > >
> > > "Warning: Failure at t=1.905843e002. Unable to
> meet
> > > integration tolerances without reducing the step
> size
> > > below the smallest value allowed (5.551115e017)
> at
> > > time t."
> > >
> > > Can anyone tell me what I need to do? How in
> general
> > > must the problem be formulated to return any
> number
> > > of internally calculated variables?
> > >
> > > Thanks!
> >
> > Is EMF = 0 at t=0 ?
> > Otherwise you should try to start with consistent
> > initial conditions, i.e. y(4,1)t=0 = EMFt=0.
> >
> > Another source for the difficulties encountered
> > could be that y(4,1) is not continuously
> > differentiable because of the numerical derivative
> > you take for psidot.
> >
> > A method that always works to get the EMFvalues at
> the
> > output times is to calculate EMF from the yvector
> > _after_ MATLAB has returned from ODE15s (in the
> same
> > manner as you do in your function ydot = feaode(t,
> y)).
> >
> > Best wishes
> > Torsten.
>
> Thank you Torsten, I have tried entering a more
> accurate calculated initial value for
> y(4,1) without much success (I did it in debug mode,
> so I didn't have the code handy to show), so perhaps
> the problem is caused by the derivative of psidot not
> being sufficiently smooth or continuous.
>
> I will do as you suggest and calculate the EMF values
> afterwards but I mainly wanted these values to check
> the system was behaving as I expected. I wanted to
> check off the EMF values against the y(3,1) values to
> ensure they made physical sense. I'm not sure working
> backwards will do this as well.
>
At t=tend, you have the solutions y(1,1), y(2,1) and
y(3,1) in Y2 for all times of T2.
If you then calculate EMF for all times in T2
with the corresponding solution values in Y2,
you get the same information as if you had used
EMF as a 4th solution variable.
> Thanks you for your help though, at least I know I am
> setting up the mass matrix correctly and that the
> problem lies elsewhere.
Best wishes
Torsten.
