clear all; close all; syms x
Define (and plot) the function (dependent variable) Symbolic functions (of one or several symbolic variables) can be defined and plotted as follows.
y = (0.65-0.01*x)*200*exp(0.025*x)-0.45*x; ezplot(y,[0,40]);
Differentiate the symbolic expressions
dydx = diff(y,x)
dydx = -2*exp(1/40*x)+1/40*(130-2*x)*exp(1/40*x)-9/20
Find the critical points of y=y(x). The solve command attempts to find explicit expressions for the roots of the equation dydx=0. This is in many cases impossible. If symbolic solutions cannot be determined, numeric solutions are returned. (See below)
xmax = solve(dydx); xmax = double(xmax)
xmax = 19.47 -107.61
Notice that we get two critical points (roots of dydx=0). We must pick only values that have a meaning in our original problem, in this case the positive value only.
xmax =xmax(1)
xmax = 19.47
Compute the numeric value ymax=y(xmax)
ymax=subs(y,x,xmax)
ymax = 139.39
Alternately, numeric solutions of dydx=0 can be obtained through Newton's Method:
F = dydx; F1 = diff(F,x); format long N = 10; % number of iterations x0 = 19 % initial guess fprintf(' iteration xvalue\n\n'); for i=1:N x1=x0-subs(F,x,x0)/subs(F1,x,x0); fprintf('%5.0f %1.16f\n', i, x1); x0 = x1; end
x0 = 19 iteration xvalue 1 19.4741582256645000 2 19.4681604158219840 3 19.4681594441612480 4 19.4681594441612300 5 19.4681594441612200 6 19.4681594441612270 7 19.4681594441612160 8 19.4681594441612230 9 19.4681594441612200 10 19.4681594441612270
display('Hence, the critical point (solution of F=0) is (approx)'), x1
Hence, the critical point (solution of F=0) is (approx) x1 = 19.46815944416123
Sensitivity to the parameter r (p. 10)
It turns out that, in this particular case, one could AGAIN solve 'explicitly' the equation dydx=0 when r is treated as a symbolic variable. This is though plain luck, since there exists a special function exists (and is programmed in the computer) that describes the solution of dydx=0. Attempting to go explicit when studying sensitivity is highly NOT RECOMMENDED way to go!
Sensitivity to the parameter r (p. 10)
xvalues = 0; for r = 0.008:0.001:0.012 y = (0.65-r*x)*200*exp(0.025*x)-0.45*x; dydx=diff(y,x); xmaxr=solve(dydx,x); xmaxr = double(xmaxr); xmaxr = xmaxr(1); xvalues = [xvalues; xmaxr]; end xvalues = xvalues(2:end);
Display the optimal times x for several values of c around 0.025
rvalues = [0.008; 0.009; 0.01; 0.011; 0.012];
format short;
display([rvalues,xvalues])
ans = 0.0080 36.7625 0.0090 27.1497 0.0100 19.4682 0.0110 13.2103 0.0120 8.0309
Plot the xmax as a function of r, indicating the values shown in the table above
figure, plot(rvalues,xvalues, '--mo',... 'LineWidth',2,... 'MarkerEdgeColor','k',... 'MarkerFaceColor',[.49 1 .63],... 'MarkerSize',12) xlabel('r'), ylabel('xmax')
Sensitivity to the parameter c (p. 10)
xvalues = 0; for c = 0.022:0.001:0.028 y = (0.65-0.01*x)*200*exp(c*x)-0.45*x; dydx=diff(y,x); xmaxc=solve(dydx); xmaxc = double(xmaxc); xmaxc = xmaxc(1); xvalues = [xvalues; xmaxc]; end xvalues = xvalues(2:end);
Display the optimal times x for several values of c around 0.025
cvalues = 0.022:0.001:0.028; cvalues=cvalues'; % transposes the row into a column format short; display([cvalues,xvalues])
ans = 0.0220 11.6263 0.0230 14.5159 0.0240 17.1166 0.0250 19.4682 0.0260 21.6037 0.0270 23.5507 0.0280 25.3322
Plot the xmax as a function of r, indicating the values shown in the table above
figure, plot(cvalues,xvalues, '--mo',... 'LineWidth',2,... 'MarkerEdgeColor','k',... 'MarkerFaceColor',[.49 1 .63],... 'MarkerSize',12) xlabel('c'), ylabel('xmax')