Question
Is the steady state stable or unstable? Why?
Write a function rhs(y::Vector{Number}) that takes in the vector $\boldsymbol{y}$ and returns the right-hand side of the two differential equations provided above, collected into a vector. Your function should use the values from the sliders above, i.e. K, G, H, U, and T_inf.
Now that we have determined the steady state solution, let us investigate the stability numerically. Call the jacobian function defined above, using either your fixed-point or Newton–Raphson solution.
Using the functions defined below, which correctly implement fixed-point iteration and the Newton–Raphson method, solve for y_ss.
Check Your Work!
If your rhs function is implemented properly, then using the default values provided above you should find rhs([0.5, 0.8]) = [-0.161873, 0.362375], as checked below
Determine the steady-state value, according to the fixed-point method, and save it in a variable titled y_ss_fp.
The function fixed_pt(f, y0) defined below solves for the vector y where f(y) = y. However, we are interested in determining where rhs(y) = 0. Thus, we need to define a new function that can be passed into the fixed_pt function. Recalling your written problem set from last week, write a function fp_rhs(y) for which y = fp_rhs(y) is a solution to rhs(y) = 0. There are many options to choose from, but not all of them will lead to a converged solution.
We are now interested in solving for the steady-states of the highly nonlinear functions above, using numerical methods. In particular, we would like to find the vector y_ss such that rhs(y_ss) = 0.
Now determine the same steady-state solution using the function newton_raphson defined below. This function takes a function f and an initial guess y0, and solves for y such that f(y) = 0.
Here, we will continue analyzing the scenario described in the written portion of the problem set. Recall that in terms of the dimensionless concentration and temperature, we have the coupled equations
$$\dfrac{\text{d} c^*}{\text{d} t^*} \, = \, 1 \, - \, c^* \, - \, \tilde{K} \, (c^*)^2 \, e^{-\tilde{G}/T^*}$$
and
$$\dfrac{\text{d} T^*}{\text{d} t^*} \, = \, 1 \, + \, \tilde{U} \, T^*_\infty \, - \, (1 + \tilde{U}) \, T^* \, + \, \tilde{K} \, \tilde{H} \, (c^*)^2 \, e^{-\tilde{G}/T^*} ~.$$
In order to be able to solve these equations numerically, we want to express them in a form convenient to a computer. To this end, we define
$$y_1 \, \equiv \, c^* \qquad \text{and} \qquad y_2 \, \equiv \, T^* ~,$$
such that our collective unknowns are contained in the vector
$$\boldsymbol{y}(t) \, = \, \begin{bmatrix} y_1 (t) \\[3pt] y_2 (t) \end{bmatrix} ~.$$
Note that in our code, the "time" variable is understood to be the dimensionless time $t^*$.
First, let us define the relevant constants in the problem. We will drop the "$\sim$" accent for notational simplicity. You can use the sliders to change the values of the different constants. The slider value is printed for your convenience.
[Type your answer here]