Fourth-Order Autonomous Ordinary Differential Equation(ODE) with Inputs (version 1.4)

This page/program provides numerical solution of a 4th-order ODE in the form of 4 equations:
          dy0/dt=f0(y0,y1,y2,y3,t)=F0(y0,y1,y2,y3)+K0*u(t),
          dy1/dt=f1(y0,y1,y2,y3,t)=F1(y0,y1,y2,y3)+K1*u(t),
          dy2/dt=f2(y0,y1,y2,y3,t)=F2(y0,y1,y2,y3)+K2*u(t),
          dy3/dt=f3(y0,y1,y2,y3,t)=F3(y0,y1,y2,y3)+K3*u(t),
          at t=0, y0(t)=y0(0), y1(t)=y1(0), y2(t)=y2(0), y3(t)=y3(0).

In math equation style:
${\displaystyle \frac{{dy}_0}{dt}=f_0(y_0,y_1,y_2,y_3,t)=F_0(y_0,y_1,y_2,y_3)+K_0\cdot u\left(t\right),}$
${\displaystyle \frac{{dy}_1}{dt}=f_1(y_0,y_1,y_2,y_3,t)=F_1(y_0,y_1,y_2,y_3)+K_1\cdot u\left(t\right),}$
${\displaystyle \frac{{dy}_2}{dt}=f_2(y_0,y_1,y_2,y_3,t)=F_2(y_0,y_1,y_2,y_3)+K_2\cdot u\left(t\right),}$
${\displaystyle \frac{{dy}_3}{dt}=f_3(y_0,y_1,y_2,y_3,t)=F_3(y_0,y_1,y_2,y_3)+K_3\cdot u\left(t\right),}$
at ${\displaystyle t=0,y_0\left(t\right)=y_0\left(0\right),y_1\left(t\right)=y_1\left(0\right),y_2\left(t\right)=y_2\left(0\right),y_3\left(t\right)=y_3\left(0\right)}$.

Another form of 4th-order ODE, ${\displaystyle y^{\left(4\right)}=F(y,y\prime ,y\prime \prime ,y\prime \prime \prime )+K\cdot u\left(t\right)}$, is a special case of the above 4-equation form. 
Let ${\displaystyle y=y_0,y\prime =y_1,y\prime \prime =y_2,,y\prime \prime \prime =y_3,K_0=K_1=K_2=0,K_3=K,}$ 
then ${\displaystyle F_0(y_0,y_1,y_2,y_3)=y_1,F_1(y_0,y_1,y_2,y_3)=y_2,F_2(y_0,y_1,y_2,y_3)=y_3,}$
and ${\displaystyle F_3(y_0,y_1,y_2,y_3)=F(y,y\prime ,y\prime \prime ,y\prime \prime \prime )=F(y_0,y_1,y_2,y_3)}$

When you open this page, the solution/simulation of a default/pre-defined example ODE is ready to run.
Click "Solve" button to display the solution for the ODE over a time period: [0, T], where T=N*dT.
Or Click "Start/Reset" button to run the ODE simulation continuously, you can stop, continue, or reset the simulation process.
A set of ODE examples are provided, use the selection menu to select pre-defined functions for the example ODE, and then stop the running ODE and start/reset again to run the new example.
Click on the "Phase Plane" equal to: (1)-Set the respective initial values as the clicked point position. and (2)-Same as click "Solve" button to start a new ODE solution process with the new initial values.

The default ODE example represents a 2-D satellite motion in gravitational force field of an imaginary planet. In  celestial mechanics, this example is a simplified model for the so-called two-body Kepler problem. The solution of this ODE is called a Kepler orbit, which can be an ellipse, a parabola, an hyperbola, or even a straight line. As an ODE solution, only an ellipse orbit can be stable and periodic, determined by a set of initial position and velocity combinations.
The variables (y0, y2) represent position, (y1,y3) represent velocity. For the phase plane plotter, select y0 and y2 to display the satellite position trajectory, and select y1 and y3 to display the velocity trajectory.
The trajectory shape and stability of the satellite are affected by the initial velocity and position. You can try different combinations of initial conditions and N and dT parameters to see their effects on the plotted ODE solutions.
Example 8 represents two coupled van der Pol oscillators. Other examples are similar to those in the ODE3 page.

To define your own ODE system, just input or re-define the ODE functions and parameters:
Functions: F0(y0,y1,y2,y3), F1(y0,y1,y2,y3), F2(y0,y1,y2,y3), and F3(y0,y1,y2,y3).
Initial values: y0(0), y1(0), y2(0), and y3(0).
K0, K1, K2, and K3, to control the input signal u(t) to each ODE channel/equation.
dT--Integration/sampling time step-size, N--Number of data to display in the 2 plotters.
Then, stop the running ODE, click "Solve" or "Start/Reset" button again to start your own ODE solution/simulation.
A more detailed usage of this ODE tool is at the end of this page.

For definition of the functions Fi(y0,y1,y2,y3,t)(i=0,1,2,3), see "Supported Operators and Functions".

Push "Start" or "Solve" Button to Run the Predefined Example

More Usage Information for the ODE Tool:

Two types of solutions are provided:
Numerical solution over a definite time period [0, T], where T=N*dT.
Run calculation/simulation/integration of the ODE continuously, to display:
the system outputs which are excited by periodical input signals,
or outputs which form limit-cycles of the ODE,
or outputs which are chaotic solutions or attractors of the ODE.
The classic explicit fourth-order Runge-Kutta method(RK4) is used as numerical integration method in this program. 

Usually the value of N can be set between 1000 to 20000. A bigger N will display more data in the graph plotter, but will also cost more memory or even cause a stack overflow problem, and hence stop the program.
Currently, the ODE solver/tool employs a 'fixed step size' integration method and the step size dT must be set by the user. Theoretically, a proper value of step size dT can be determined according to the ODE's maximum state changing rate, maximum inherent frequency(or minimum period), minimum time-constant, and/or the so-called stiffness.
As a practical approach, you can begin with a small dT and use the "Solve" button, to get initial yi(t)(i=0,1,2,3) displays, and then use the yi(t) displays to obtain a rough estimate of the ODE's minimum time-constant(or minimum period).
Usually dT can be set between 1/20 to 1/5 of the minimum time-constant (or minimum period) for most non-stiff ODEs. A smaller dT will lead to higher accuracy, but will cause longer computation time as well.
You can try different combinations of dT and N values, to get a stable and accurate ODE solution, and also nice graph displays, at the same time.
A special tool for exploring relationship between the step size and accuracy of the ODE solution will be developed and added to our website.
In the future, the so-called 'adaptive step size' algorithms may be added to our ODE solvers/tools. Then, the step size will be variable during the integration process according to the dynamical characteristics of the ODE. That is, the step size dT will be automatically determined by the ODE solver itself, and the user just need to set error tolerance parameters optionally.

The input signal u(t) is provided as a set of selectable time signals for test and verification purposes.
Currently, eight types of input signals u(t) are provided: five types of periodic signals, one constant signal, and two types of stochastic signals.
For periodic u(t), you can set signal period Ts, and magnitude Am. To keep accuracy of the periodic u(t), the value of Ts is limited to Ts>=50*dT.
Set amplitude Am=0 will let u(t)=0, same as remove the input signal. Ki(i=0,1,2,3) are used to apply u(t) to channel i. Set Ki=0 will remove the input signal from the respective channel, and Ki>1 will amplify the signal.

The user can select 2 of 4 yi(t)(i=0,1,2,3) variables for phase plane plotter. After phase plane coordinates are established, move the Mouse cursor on the phase plane will display the 2 yi(t) position. Click the Mouse left key, will start a "Solve" process with the displayed 2 yi(t) as the initial condition of the ODE, that is, the solution trajectory will start from the yi(t) point. Notice that both the phase plane and oscilloscope plotters are auto-scaled and auto-centered, the clicked initial point may be displayed in a new coordinate grid immediately. Before click again, you need to move the Mouse to update the mouse position display on the phase plane.
The phase plane graph is auto-scaled and auto-centered on the canvas. If the checkbox "Keep Graph Shape or Not" is checked, the plotter uses the same scale factor in both x and y directions and the graph will keep its original shape. If the checkbox is unchecked, the graph will be stretched in one direction to make full use of the canvas area.
On the phase plane graph, a series of colored dots, from blue to red, are used to indicate the moving direction and/or moving speed of the solution points on the plane.
If you changed any input parameter or the function selection menu, you need to "stop" the program and "start" again to let the new parameter and/or function selection be read and used in the solution.
For "Solve" button started solution, if you changed any checkbox status, you also need to click "Solve" button again to start a new solution with the new checkbox status.  
In case an error occurs, the program will display a notice and may indicate the cause of the error. After the error is corrected, Click "Solve" or "Start/Reset" button to re-start a new solution process.