Euler with step hand written and Matlab solution

Question

A biological process includes the growth of biomass (B) from substrate (S). The material balances on this batch process yield:

where B and S are the respective biomass and substrate concentrations. The reaction kinetics are described by kr= 0.70000 and K= 0.00300 in consistent units. At the beginning (t=0) the initial values in consistent units are (B,S)=(0.05, 5).

1. Use Euler's integration performing two steps with dt=1 , to find the values of B and S at t=1 and t=2. Describe in detail all calculations by hand.
2. Use one step of 4th order Runge-Kutta with dt=1 to find the values of B and S at t=1. Describe in detail all calculations by hand and compare to your findings with Euler's.
3. Use MATLAB and ode45 to extract the dynamic behavior up to t=50. Provide a plot with B and S as a function of time. Do the results appear reasonable?
4. Use MATLAB and ode15s to extract the dynamic behavior up to t=50. Provide a plot with B and S as a function of time. Do the results appear reasonable? Contrast to the findings in the previous step. [ 1pts]
5. Based on the results with ode15s:
6. Provide estimates for B and S at time t=5.0 
7. Provide estimates for the time t and the value of B when S=0.001
8. As part of stability analysis, provide the array of functions describing the Jacobian of the system.
9. Find the eigenvalues of the Jacobian at t=5.0 and the second time-point where S=0.001. Contrast the values of the eigenvalues at these two points in time.

Answer 1

Hi,

Because this question looks like an exact question for an assignment, I'll only partially answer the question, skipping the hand calculations and focusing primarily on MATLAB/Octave code. I'll also only focus on 2 of the integration methods: Runge-Kutta (RK) and ode15s - again in keeping with academic honesty requirements. The parameter N (number of steps) can be changed to 2 to answer the second question.

In the script below, both methods are implemented. It is clear that the system is unstable at high values of N. A 4th order RK integration is carried out in a for loop. For ode15s, the option settings ("opt = " line) are required to stabilize the integration algorithm. Without these settings, the algorithm ends up in a NULL situation. A relative and absolute tolerance of 1e-2 is pretty relaxed, so it is a good idea to always put at least this high a value in (make it smaller for higher accuracy).

Unfortunately, the platform doesn't give me a way to post the resulting image. However, if you run the script below, you'll see that the B parameters increases exponentially and the S parameter decreases in a similar way. From the first and second equation, its also pretty clear that dB_dS is a larger negative constant number.

BR,

Omar

----------------------------------------------------------

clear;

close all;

% t0, S0, and B0 are the initial conditions

t0 = 0; B0 = 0.05; S0 = 5;

% kr and K are the system parameters

kr = 0.7; K = 0.003;

% system is defined by the following equations

% MATLAB allows you to define functions as dependent equations

% in this way.

dB_dt = @(B,S) kr*B*S/(K+S);

dS_dt = @(B,S) -0.75*kr*B*S/(K+S);

f = @(t,y) [dB_dt(y(1),y(2)); dS_dt(y(1),y(2))];

% visualization

% I find it helps to use these commands to set up the figure

figure(1); hold on; grid on;

% ------ Runge-Kutta integration -----------

% This integration is broken down into 4 parts (for 4th order)

% N is the number of integration steps

% dt is the step size

dt=1; N=50;

% these are the initial conditions

t(1) = t0; y(:,1) = [B0; S0];

for i = 1:N

% the 4 components of 4th order RK integration

k1 = f(t(i), y(:,i));

k2 = f(t(i)+dt/2,y(:,i)+dt/2*k1);

k3 = f(t(i)+dt/2,y(:,i)+dt/2*k2);

k4 = f(t(i)+dt,y(:,i)+dt*k3);

y(:,i+1) = y(:,i) + dt*(k1+2*k2+2*k3+k4)/6.0;

t(i+1) = t(i) + dt;

end

% this is to display the results. it might be helpful to

% reduce N to 10 to take a good look at the output screen.

disp(' --- Runge-Kutta --- ');

disp([t; y]); % print values in the command window

plot(t, y); % plot B and S (now contained in y) with respect to time t

% ------ ode15s integration -----------

y0 = [B0; S0];

tol = 1e-2;

opt = odeset('RelTol',tol, 'AbsTol',tol);

[t,y] = ode15s(f,[0, dt*N], y0, opt);

disp(' --- ode15s --- ');

disp([t'; y']); % print values in the command window

plot(t, y); % plot B and S (now contained in y) with respect to time t

ylim([-1000 1000]);

% finalize plot outputs

% set the x and y labels and the legend

xlabel('Time, s');

ylabel('Process Value, -');

legend('R-K B', 'R-K S', 'ode15s B', 'ode15s S');