Note that there are some explanatory texts on larger screens.

plurals
  1. PO
    primarykey
    data
    text
    <p>You said that you need to find the root "with great accuracy." The answer by @CST-Link above is not how to do this if you're solving the differential equation numerically with <code>ODE45</code>. In fact it's a bad idea. It requires outputting your solution points at high resolution and introduces error, which can be particularly bad if the equations you are integrating are conservative (i.e., you will effectively add or subtract energy from a solution where the total energy should always remain constant). You need to use the event detection (zero crossing) functionality of Matlab's ODE solvers. This will generally be able to find your root with an accuracy close to machine epsilon, <code>eps</code>. Whereas the technique given by @CST-Link will only give an accuracy on the order of ODE45's step-size, which can be very large (interpolation can be used to help improve this, but you're not going to get close to <code>eps</code> unless you use tiny step sizes).</p> <p>Looking at Matlab's help for <code>ODE45</code>, the events functionality can seem confusing, so I'll try to give a simpler example using code based on <code>ballode</code> (<code>help ballode</code> and <code>edit ballode</code> for more info), an example that Matlab includes to demonstrate events.</p> <pre><code>function eventsdemo options = odeset('Events',@efun); % specify name of events function [t,y,te,ye,ie] = ode45(@f,[0 10],[0;20],options); % integrate figure plot(t,y(:,1),'b',te,ye(:,1),'r.') function dydt = f(t,y) dydt = [y(2);-9.8]; % differential equation for ballistic motion function [value,isterminal,direction] = efun(t,y) value = y(1); isterminal = 1; direction = -1; </code></pre> <p>The example is just simple ballistic motion of a ball in the vertical direction: the function <code>f</code>. The goal is to accurately detect the velocity of the ball at the time it crosses the ground plane, <code>y(1) == 0</code> line in the negative direction. If you detect contact too soon the ball's velocity will be less than reality and energy will be lost; too late and energy will be injected. The <code>value</code> output of <code>efun</code> defines the equation that must equal zero for the crossing. This can be as simple or as complex as needed, can depend on all of the state variables and time, and you can detect multiple crossings by specifying a vector. In your case, it sounds like you are interested in roots along the x-axis, so I imagine that you might have an identical definition of <code>value</code>. If you only want the first root or if there's only one, then <code>isterminal</code> will stop the integration when it is found. Finally, if you don't know the slope/gradient of your function at the root, you can set <code>direction</code> to zero. Try out this event function with the above code (<code>ie</code> can be used to tell which of the three events is triggered in the <code>te</code> and <code>ye</code> outputs):</p> <pre><code>function [value,isterminal,direction] = efun(t,y) value = [y(2) y(1)-10 y(1)]; isterminal = [0 0 1]; direction = [-1 0 -1]; </code></pre> <p>Some minor caveats. You need to set the final integration time to be sufficiently long. In other words, the event must occur between <code>t0</code> and <code>tf</code> (see <code>ballode</code> for a scheme to iterate calls to <code>ODE45</code> if you don't know where in time your events will occur). If you set <code>isterminal</code> to zero, the output <code>t</code> vector and <code>y</code> matrix will be trimmed to end at the time of the terminal event. Lastly, if you specify fixed step-size outputs, e.g., <code>TSPAN = t0:dt:tf</code>, the last output point will likely have a smaller step-size.</p>
    singulars
    1. This table or related slice is empty.
    plurals
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. COThis is a very informative solution, and it taught me a great deal about the options of `ode45`, stuff that I didn't know and for which I thank you. But, in defense of my post, :-) I'd like to state again the final question that was asked (and, with your permission, I will quote): *So, does anyone know how I could find a "root" to a set of given Y values without knowing the actual equation?*
      singulars
      1. This table or related slice is empty.
    2. COI assume that you at least have the differential equations and that you simulated them yourself? I interpreted your question as you not having the an analytical solution to the differential equations (what you call "the actual equation" of the Y values), hence the reason for using `ode45` to obtain a numeric solution. If this is not true, you need to rewrite your question and remove reference to `ode45` - the means for obtaining the data would be irrelevant. My solution is the proper way of numerically finding roots of differential equations where an analytical solution can't be found.
      singulars
    3. COSpeaking of ODE solver options, if you use the single output argument form with your numerical integration function you get a struct that can be used with the `deval` function after the the fact to interpolate the solution. The interpolation is "optimal" in that the [interpolating polynomial](http://en.wikipedia.org/wiki/Polynomial_interpolation) of the ODE solver is used (`edit private/ntrp45`). This won't necessarily give as good results as with events, but it is the best option if you need finer-grained output after applying some form of @CST-Link's solution.
      singulars
 

Querying!

 
Guidance

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload