Math 448/548 - Example 3.1. Nonlinear Pig (p.61-71)

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')