# Simulating the response of nonlinear system with stiff differential equations

I want to simulate the response of a nonlinear system given in the following form:

$$\dot{x_1} = f_1(\bar{x_1})+g_1(\bar{x_1})x_2, \ x_1(0) = 0.2$$ $$\dot{x_2} = f_2(\bar{x_2})+g_2(\bar{x_2})x_3, \ x_2(0) = 0.5$$ $$\dot{x_2} = f_3(\bar{x_3})+g_3(\bar{x_3})u, \ x_3(0) = 0.5$$ $$y = x_1$$

where $$u$$ is the control input to the system, $$y$$ is the outputof the systen and the functions $$f_i$$ and $$g_i$$, $$i = 1,2,3$$ are known and nonlinear. In order to come up with the control input at each time step some computations have to be carried out, which are the following:

$$\epsilon_1 = S_1^{-1}(\frac{x(1)-y_d}{p_1(t)})$$ $$n_1 = k_1[\frac{1}{p_1(t)\frac{dS_1}{d\epsilon_1}(\epsilon_1)}+p_1(t)\frac{dS_1}{d\epsilon_1}(\epsilon_1)]\epsilon_1$$ $$\dot{ζ_1} = \frac{\epsilon_1n_1}{p_1(t)\frac{dS_1}{d\epsilon_1}(\epsilon_1)}$$ $$a_1 = N(ζ_1)n_1$$

$$\epsilon_2 = S_2^{-1}(\frac{x(2)-a_1}{p_2(t)})$$ $$n_2 = k_2[\frac{1}{p_2(t)\frac{dS_2}{d\epsilon_2}(\epsilon_2)}+p_2(t)\frac{dS_2}{d\epsilon_2}(\epsilon_2)]\epsilon_2$$ $$\dot{ζ_2} = \frac{\epsilon_2n_2}{p_2(t)\frac{dS_2}{d\epsilon_2}(\epsilon_2)}$$ $$a_2 = N(ζ_2)n_2$$

$$\epsilon_3 = S_3^{-1}(\frac{x(3)-a_2}{p_3(t)})$$ $$n_3 = k_3[\frac{1}{p_3(t)\frac{dS_3}{d\epsilon_3}(\epsilon_3)}+p_3(t)\frac{dS_3}{d\epsilon_3}(\epsilon_3)]\epsilon_3$$ $$\dot{ζ_3} = \frac{\epsilon_3n_3}{p_3(t)\frac{dS_3}{d\epsilon_3}(\epsilon_3)}$$ $$u = N(ζ_3)n_3$$

The functions $$S_i(\epsilon_i)$$, $$i=1,2,3$$ are some smooth bijective functions, which are known and so their inverse functions and their derivatives can be computed. Same goes for the functions $$p_i(t)$$, $$i=1,2,3$$ which are some performance functions. Finally, the function $$N(ζ)$$ is an even Nussbaum function, also known and $$y_d(t)$$ is a desired trajectory for the output of the system. The numerical solution of the closed loop (with $$u$$ as the input) is stiff and so I will use the ode15s MATLAB solver for the simulation. Now, my problem is how can I compute the values of the $$ζ_i$$, $$i=1,2,3$$ variables in order to compute the value of $$u$$. Some pseudocode I have thought off is:

function dx = odefun(t,x)

compute ε1, p1, n1, ζ1
compute a1

compute ε2, p2, n2, ζ2
compute a2

compute ε3, p3, n3, ζ3
compute u

dx = [x1_dot ; x2_dot ; x3_dot]

end


I know I can compute the values of the variables $$\epsilon_i, n_i, p_i$$ by using some function handles but know comes my problem regarding the variables $$ζ_i$$ due to the fact that they are computed though the differential equations. Is there a way to include them in the ode solver and solve them all at once in the proper way ?

• Is there a reason you can't just consider this as a system of 6 coupled ODEs? Jun 4 '20 at 16:53
• @whpowell96 Because in order to compute the intermediate control signals $a_i$, I need to compute first the values of $ζ$ and then continue. I don't know if I can do it all in one and that's why I am asking. Jun 4 '20 at 17:13
• If you consider this as 6 ODEs for the variables $x_1,x_2,x_3,\zeta_1,\zeta_2,\zeta_3$, then you will always "know" the values of $\zeta$at the same timestep as $x$ since they are computed simultaeneously by the ODE solver Jun 4 '20 at 19:38