University of Florida/Egm6341/s10.Team2/HW7

Statement

Pg.38-4 Show the case of stable growth for Verhulst equation

Solution

Author

Egm6341.s10.team2.lee 20:56, 23 April 2010 (UTC)

Statement

P.39-1

Solution

We have Verhulst model or the Logistic equation as P.38-3

     ${\displaystyle {\dot {x}}={\frac {dx}{dt}}=rx(1-{\frac {x}{x_{max}}})}$

Separating variables and integrating we have

     ${\displaystyle \int _{x_{0}}^{x}{\frac {x_{max}}{x(x_{max}-x)}}dx=\int _{t_{0}=0}^{t}rdt}$

which can be written as

     ${\displaystyle \int _{x_{0}}^{x}{\frac {1}{x}}dx+\int _{x_{0}}^{x}{\frac {1}{(x_{max}-x)}}dx=\int _{t_{0}=0}^{t}rdt}$

We get

     ${\displaystyle log_{e}({\frac {x}{x_{0}}})-log_{e}({\frac {x_{max}-x}{x_{max}-x_{0}}})=r(t-0)}$

     ${\displaystyle {\frac {x(x_{max}-x_{0})}{x_{0}(x_{max}-x)}}=e^{rt}}$

Rearranging we have

     ${\displaystyle x={\frac {x_{max}x_{0}e^{rt}}{x_{max}+x_{0}(e^{rt}-1)}}}$

Author

--Egm6341.s10.team2.niki 12:31, 23 April 2010 (UTC)

Statement

Pg.39-1

Use the Hermit-Simpson Algorithm to integrate the Verhulst population growth model on Pg.38-3:

${\displaystyle {\frac {\mathrm {d} x}{\mathrm {d} t}}=rx\left(1-{\frac {x}{x_{max}}}\right)}$

Consider the following: ${\displaystyle x_{max}=10\;\;r=1.2\;\;t\in \left[0,10\right]}$ and
Case 1: ${\displaystyle \;x_{0}=2}$
Case 2: ${\displaystyle \;x_{0}=7}$

And also plot the two cases using matlab.

Solution

Matlab code

%Matlab code to integrate the the Verhulst population growth model
clear;
clc;
xmax=10;
r=1.2;
t=0:0.1:10;
h=0.1; %time step size
syms x;
syms x11;
f=r*x*(1-(x/xmax));
%case 1: when x0=2
for i=1:1:100
    x1(1)=2;
    f1=subs(f,x,x1(i));
    f11=subs(f,x,x11);
    x1half=((1/2)*(x1+x11))+((h/8)*(f1-f11));
    f1half=subs(f,x,x1half);
    F=x1(i)-x11+(h/6)*(f1+f11+(4*f1half));
    guessx11(i)=x1(i);
    dummy=1;
    while dummy==1
    newx11(i)=guessx11(i)-((1/(subs(diff(F(i)),x11,guessx11(i))))*(subs(F(i),x11,guessx11(i))));
    if double(abs(newx11(i)-guessx11(i)))<=10^(-6)
       x1(i+1)=newx11(i);
       dummy=0;
        break;
    else
        guessx11(i)=newx11(i);
        dummy=1;
    end
    end
end
plot(t,x1);

THE PLOT SHOWING THE POPULATION CHANGE W.R.T TIME, HERE VERHULST MODEL IS USED TO COMPUTE THE POPULATION WHERE XO=INITIAL POPULATION=2

THE PLOT SHOWING THE POPULATION CHANGE W.R.T TIME, HERE VERHULST MODEL IS USED TO COMPUTE THE POPULATION WHERE XO=INITIAL POPULATION=7

Author

Solved by--Srikanth Madala (SM)

Statement

Pg.40-1

Produce the diagrams showing the bifurcation diagram for the logistic map and its sensitivity to initial conditions, as shown on figures 15.6 and 15.7, of the book "Differential Equations: Linear,Nonlinear,Ordinary,Partial" by King,Billingham and Otto (Pg.455-456).

Solution

for r= 0:0.005:4
    x=rand(1);
    for j = 1:100
        x=r*x*(1-x);
    end
    xout=[];
    for j=1:400
        x=r*x*(1-x); xout=[xout x];
    end
    plot (r*ones(size(xout)),xout,'.','MarkerSize',3)
    axis([0 4 0 1]), hold on, pause(0.01)
end

function [xn1 xn2]=fig157

r=4;

format long

x01=.1;
x02=0.1+(10^-16);

for j=0:100
k=j+1;
x1=r*x01.*(1-x01);
x2=r*x02.*(1-x02);

xn1(k,1)=x1;
xn2(k,1)=x2;

x01=x1;
x02=x2;

end

n=0:100;
plot(n,xn1,'o')
hold on
plot(n,xn2,'x')

Author

--Egm6341.s10.Team2.GV 16:28, 23 April 2010 (UTC)

Statement

Pg.40-2

Solution

Matlab code

clear;
clc;
xmax=10;
r=1.2;
t=0:0.1:10;
h=0.1; %time step size
syms x;
syms x11;
f=r*x*(1-(x/xmax));
%case 1: when x0=2
for k=1:1:5
    h=2*h;
for i=1:1:100
    x1(1)=2;
    f1=subs(f,x,x1(i));
    f11=subs(f,x,x11);
    x1half=((1/2)*(x1+x11))+((h/8)*(f1-f11));
    f1half=subs(f,x,x1half);
    F=x1(i)-x11+(h/6)*(f1+f11+(4*f1half));
    guessx11(i)=x1(i);
    dummy=1;
    while dummy==1
    newx11(i)=guessx11(i)-((1/(subs(diff(F(i)),x11,guessx11(i))))*(subs(F(i),x11,guessx11(i))));
    if double(abs(newx11(i)-guessx11(i)))<=10^(-6)
       x1(i+1)=newx11(i);
       dummy=0;
        break;
    else
        guessx11(i)=newx11(i);
        dummy=1;
    end
    end
end
plot(t,x1);
hold on
end

THE PLOT SHOWING THE POPULATION CHANGE W.R.T TIME, HERE VERHULST MODEL IS USED TO COMPUTE THE POPULATION WHERE XO=INITIAL POPULATION=2

Author

Solved by--Srikanth Madala (SM)

Statement

Pg.40-3 and Pg.41-1

Solve the following Group of ODE's as they relate to the problem mentioned in the lecture pages mentioned above. The problem is of the motion of an airplane through different maneuvers.

for a vector with the following parameters:

${\displaystyle z=[x,y,V,\gamma ]}$
${\displaystyle {\dot {z}}=[{\dot {x}},{\dot {y}},{\dot {V}},{\dot {gamma}}]}$

${\displaystyle {\dot {x}}=Vcos(\gamma )}$

${\displaystyle {\dot {y}}=Vsin(\gamma )}$

${\displaystyle {\dot {v}}=({\frac {T-D}{m}}cos(\alpha )+{\frac {L}{m}}sin(\alpha )-gsin(\gamma )}$

${\displaystyle {\dot {\gamma }}=({\frac {T-D}{mV}}sin(\alpha )+{\frac {L}{mV}}cos(\alpha )-{\frac {g}{V}}cos(\gamma )}$

Initial Conditions are as follows:

${\displaystyle \gamma _{0}=0\quad degrees}$

${\displaystyle V_{0}=272\quad meters/sec}$

${\displaystyle x_{0}=0\quad meters}$

${\displaystyle y_{0}=0\quad meters}$

The Physical Modeling Parameters to be used are

${\displaystyle m=1005\quad Kg}$
${\displaystyle g=9.81\quad m/sec^{2}}$
${\displaystyle Sref=0.3376\quad m^{2}}$
${\displaystyle A1=-1.9431\qquad A2=-0.1499\qquad A3=0.2359}$
${\displaystyle B1=21.9\qquad B2=0}$
${\displaystyle C1=3.312x10^{-9}\quad Kgm^{-5}}$
${\displaystyle C2=-1.142x10^{-4}\quad Kgm^{-4}}$
${\displaystyle C3=1.1224\quad Kgm^{-3}}$

${\displaystyle \rho =C1y^{2}+C2y+C3}$
${\displaystyle L={\frac {1}{2}}C_{l}\rho V^{2}Sref}$
${\displaystyle C_{l}=B1\alpha +B2}$
${\displaystyle D={\frac {1}{2}}C_{d}\rho V^{2}Sref}$
${\displaystyle C_{d}=A1\alpha ^{2}+A2\alpha +A3}$

The Inputs are defined as follows:

Input T:

for ${\displaystyle 0\leq t<27seconds=6000N}$
for ${\displaystyle 27\leq t<33seconds=1000N}$
for ${\displaystyle 33\leq t\leq 40seconds=6000N}$

Input ${\displaystyle \alpha }$:

for ${\displaystyle 0\leq t<21seconds=0.03}$
for ${\displaystyle 21\leq t<27seconds=0.09to0.13}$ as it varies linearly
for ${\displaystyle 27\leq t<33seconds=-0.13to-0.2}$ as it varies linearly
for ${\displaystyle 33\leq t\leq 40seconds=-0.2to-0.13}$ as it varies linearly

Solution

% Function that calculates the thrust for any given time, t
function T = thrust(t)

%for t=0:40
 %   k=t+1;
if (t >= 0) & (t < 27),
   T = 6000;
elseif (t >= 27) & (t < 33),
   T = 1000;
else
   T = 6000;
end;
%end

% Function calculates the angle of attack from 0-40secs
function alpha = angleatak(t)

%for t=0:40
 %   k=t+1;
if (t >= 0) & (t < 21),
   alpha = 0.03;
elseif (t >= 21) & (t < 27),
   alpha = ((0.13-0.09)/(27-21))*(t-21) + 0.09;
elseif (t >= 27) & (t < 33),
   alpha = ((-0.2 + 0.13)/(33-27))*(t-27) - 0.13;
else
   alpha = ((-0.13 +0.2)/(40-33))*(t-33) - 0.2;
end;
%end

    ${\displaystyle Z_{i+1}=z_{i}+{\frac {h/2}{3}}\left[f_{i}+4f(g(z_{i},z_{i+1}))+f(z_{i+1})\right]\qquad \qquad EQUATION1}$

    ${\displaystyle g(z_{i},z_{i+1})={\frac {1}{2}}(z_{i}+z_{i+1})+{\frac {h}{8}}\left[(f(z_{i})-f(z_{i+1})\right]\qquad \qquad EQUATION2}$

    ${\displaystyle F(Z_{i+1})=0=Z_{i+1}-z_{i}-{\frac {h/2}{3}}\left[f_{i}+4f(g(z_{i},z_{i+1}))+f(z_{i+1})\right]\qquad \qquad EQUATION3}$

function f = fcalc(t,z)

% z(1) = x : horizontal position
% z(2) = y : vertical position(height)
% z(3) = v : velocity
% z(4) = gamma : angle to new x-axis

% Set parameters
m = 1005; % kg
g = 9.81; % m/s^2
Sr = 0.3376; % m^2
A1 = -1.9431;
A2 = -0.1499;
A3 = 0.2359;
B1 = 21.9;
B2 = 0;
C1 = 3.312e-9; % kg/m^5
C2 = -1.142e-4; % kg/m^4
C3 = 1.224; % kg/m^3

alpha = angleatak(t); % Call alpha
T = thrust(t); % Call thrust

Cd = A1.*alpha.^2 + A2.*alpha + A3;
Cl = B1.*alpha + B2;
rho = C1*z(4)^2 +C2*z(4) + C3; % density

D = 0.5*Cd*rho*z(3)^2*Sr; % Calculate D
L = 0.5*Cl*rho*z(3)^2*Sr; % Calculate L

f = zeros(4,1);
f(1) = z(2)*cos(z(4));
f(2) = z(2)*sin(z(4));
f(3) = (T-D)*cos(alpha)/m - L*sin(alpha)/m -g*sin(z(4));
f(4) = (T-D)*sin(alpha)/(m*z(3)) + L*cos(alpha)/(m*z(3))-(g*cos(z(4)))/z(3);

function df = dfcalc(t,z)

% z(1) = x : horizontal position
% z(2) = y : vertical position(height)
% z(3) = v : velocity
% z(4) = gamma : angle to new x-axis

% Set parameters
m = 1005; % kg
g = 9.81; % m/s^2
Sr = 0.3376; % m^2
A1 = -1.9431;
A2 = -0.1499;
A3 = 0.2359;
B1 = 21.9;
B2 = 0;
C1 = 3.312e-9; % kg/m^5
C2 = -1.142e-4; % kg/m^4
C3 = 1.224; % kg/m^3

alpha = angleatak(t); % Call alpha
T = thrust(t); % Call thrust

Cd = A1.*alpha.^2 + A2.*alpha + A3;
Cl = B1.*alpha + B2;
rho = C1*z(4)^2 +C2*z(4) + C3; % density

D = 0.5*Cd*rho*z(3)^2*Sr; % Calculate D
L = 0.5*Cl*rho*z(3)^2*Sr; % Calculate L

%Caculates DF
df=zeros(4,4);
df(1,1) = 0; %
df(1,2) = 0; %
df(1,3) = cos(z(4));
df(1,4) = -z(3)*sin(z(4));

df(2,1) = 0;
df(2,2) = 0;
df(2,3) = sin(z(4));
df(2,4) = z(2)*cos(z(1));

df(3,1) = 0;
df(3,2) = -Cd*(2*C1*z(4)+C2)*z(2)^2*Sr*cos(alpha)/(2*m) - Cl*(2*C1*z(4)+C2)*z(2)^2*Sr*sin(alpha)/(2*m);
df(3,3) = -Cd*rho*z(3)*Sr*cos(alpha)/m - Cl*rho*z(3)*Sr*sin(alpha)/m;
df(3,4) = -g*cos(z(4));

df(4,1) = 0;
df(4,2) = -Cd*(2*C1*z(2)+C2)*z(3)*Sr*sin(alpha)/(2*m) + Cl*(2*C1*z(2)+C2)*z(3)*Sr*cos(alpha)/(2*m);
df(4,3) = -T*sin(alpha)/(m*z(3)^2) - Cd*rho*Sr*sin(alpha)/(2*m) + Cl*rho*Sr*cos(alpha)/(2*m) + g*cos(z(4))/z(3)^2;
df(4,4) = g*sin(z(4))/z(3);

% This function returns the solution of EOMs given
% in Subchan and Zbikowski, Optim. Control Appl. Meth. 2007; 28:311-353.
% This function uses Hermite-Simpson (HS) algorithm and Newton-Raphson (NR) method.
% Note:
% z_out(time,1) = gamma : angle b/w aircraft body and x-axis for a given time
% z_out(time,2) = v : aircraft velocity for a given time
% z_out(time,3) = x : x-coordinate of aircraft for a given time
% z_out(time,4) = y : y-coordinate of aircraft for a given time
% t_out = times corresponding to nodes
function [t,z] = NRplane

i = 0;

% This "while" loop is for increasing # of nodes to achieve the given tolerance
while (true),
   i = i + 1; % # of iterations
   n = 2^(i-1); % # of subinterval

   tinit = 0.0; % initial time (sec)
   tfinal = 40.0; % final time (sec)

   %initial conditions
   z0 = zeros(4,1);
   z0(1) = 0;
   z0(2) = 30;
   z0(3) = 272;
   z0(4) = 0;
   f0 = fcalc(tinit,z0); % Call f(z) for t = 0

   h = (tfinal - tinit)/n; % time step

   z_out = zeros(n+1,4); % output solution
   t_out = zeros(n+1,1); % output time evolution

   % assign the initial conditions to the solution for t = 0
   z_out(1,1) = z0(1); z_out(1,2) = z0(2);
   z_out(1,3) = z0(3); z_out(1,4) = z0(4);
   t_out(1) = tinit;

   % We need the values of previous time step in HS algorithm
   zj_p = zeros(4,1); % a column vector
   zj_p(1) = z0(1); zj_p(2) = z0(2);
   zj_p(3) = z0(3); zj_p(4) = z0(4);
   tj_p = tinit;

   % This "for" loop is for time evolution.
   % "j = 1" corresponds to t = tinit and "j = n+1" corresponds to t = tfinal
   for j = 2:n+1,
       tj = tinit + (j-1)*h; % current time
       zj = zj_p; % initial guess for NR method
       k = 0;
       % This "while" loop is for NR method.
       while (true),
          k = k + 1;
          fb_p = fcalc(tj_p,zj_p); % Call f(z) for t = t(j)
          fb = fcalc(tj,zj); % Call f(z) for t = t(j+1)

          tj_md = 0.5*(tj_p + tj); % t = t(i+1/2)
          g = 0.5*(zj_p + zj) + (h/8)*(fb_p - fb); % Obtain z for t = t(j+1/2)
          fb_md = fcalc(tj_md,g); % Call f(z) for t = t(j+1/2)

          % Calculate F(z) for t = t(j+1)
          ff = zj - zj_p - (h/6)*(fb_p + 4.0*fb_md + fb);

          dfb = dfcalc(tj,zj); % Call df/dz for t = t(j+1)
          dg = 0.5*eye(4,4) - (h/8)*dfb; % Calculate dg/dz for t = t(j+1)
          dfbdg = dfcalc(tj_md,g); % Call df/dz for t = t(j+1/2)

          % Calculate dF/dz for t = t(j+1)
          dff = eye(4,4) - (h/6)*(4.0*dfbdg*dg + dfb);
          inv_dff = inv(dff); % obtain inverse of dF/dz
          zj_f = zj - inv_dff*ff; % Calculte new solution z for t = t(j+1)

          % Check tolerance
          if (norm(abs(zj_f-zj),inf) < 1.e-6),
             break;
          end;

          % Redo interation if the error is still larger than tolerance
          zj = zj_f;
       end;

       % Assign the satisfied solution to the output matrices
       t_out(j) = tj;
       z_out(j,1) = zj_f(1); z_out(j,2) = zj_f(2);
       z_out(j,3) = zj_f(3); z_out(j,4) = zj_f(4);

       % Current values become previous values for the next time step
       tj_p = tj;
       zj_p = zj_f;
   end;

   % check tolerance
   if (i ~= 1) & (norm(abs(zj_f - zj_fp),inf) < 1.e-2),
      break;
   end;

   zj_fp = zj_f;
end;

Author

--Egm6341.s10.Team2.GV 23:41, 27 April 2010 (UTC)

Statement

Integrate the Logistic equation using the inconsistent Trapezoidal-Simpson approximation ${\displaystyle x_{i+1/2}=x_{i}+x_{i+1}}$

Pg.41-2

Solution

The Matlab code is a slight alteration of the code used to integrate above.

  %Matlab code to integrate the the Verhulst population growth model
  clear;
  clc;
  xmax=10;
  r=1.2;
  t=0:0.1:10;
  h=0.1; %time step size
  syms x;
  syms x11;
  f=r*x*(1-(x/xmax));
  %case 1: when x0=2
  for i=1:1:100
   x1(1)=2;
   f1=subs(f,x,x1(i));
   f11=subs(f,x,x11);
   x1half=((1/2)*(x1+x11))%Trapezoidal_simpson approxiamtion;
   f1half=subs(f,x,x1half);
   F=x1(i)-x11+(h/6)*(f1+f11+(4*f1half));
   guessx11(i)=x1(i);
   dummy=1;
   while dummy==1
   newx11(i)=guessx11(i)-((1/(subs(diff(F(i)),x11,guessx11(i))))*(subs(F(i),x11,guessx11(i))));
   if double(abs(newx11(i)-guessx11(i)))<=10^(-6)
      x1(i+1)=newx11(i);
      dummy=0;
       break;
   else
       guessx11(i)=newx11(i);
       dummy=1;
   end
   end
  end
  plot(t,x1);

for ${\displaystyle x_{0}=2}$

for ${\displaystyle x_{0}=7}$

Author

--Egm6341.s10.team2.niki 20:33, 28 April 2010 (UTC)

Statement

Pg.42-1

Solution

Author

Statement

Pg.42-1

Given the parameterization for an ellipse as follows for the x and y coordinates, Eliminate 't' as a parameter.

${\displaystyle x=a\cos(t)\qquad \qquad y=b\sin(t)}$

Solution

${\displaystyle x=a\cos(t)}$
${\displaystyle x^{2}=a^{2}\cos ^{2}(t)}$
${\displaystyle y=b\sin(t)}$
${\displaystyle y^{2}=b^{2}\sin ^{2}(t)}$

${\displaystyle {\frac {x^{2}}{a^{2}}}=\cos ^{2}(t)}$
${\displaystyle {\frac {y^{2}}{b^{2}}}=\sin ^{2}(t)}$

${\displaystyle {\frac {x^{2}}{a^{2}}}+{\frac {y^{2}}{b^{2}}}=\cos ^{2}(t)+\sin ^{2}(t)=1}$
At this point t can be eliminated as a parameter.

Author

--Egm6341.s10.Team2.GV 17:51, 23 April 2010 (UTC)

Statement

Pg.42-2

Solution

Author

Statement

Pg.42-2 To show the integral for elliptical circumference

Solution

Author

Egm6341.s10.team2.lee 20:58, 23 April 2010 (UTC)

Statement

Pg.42-2

Solution

Author

Statement

Pg.42-2

To show that ${\displaystyle I=\int _{-1}^{+1}f(x)dx=\int _{0}^{\pi }f(cos\theta )sin\theta d\theta }$

Solution

We have the integral as

${\displaystyle I=\int _{-1}^{+1}f(x)dx}$.

Using the transformation ${\displaystyle \displaystyle {x=cos\theta }}$ we have

${\displaystyle \displaystyle {f(x)=f(cos\theta )}}$

${\displaystyle \displaystyle {dx=-sin\theta d\theta }}$

and using the inverse transformation the limits ${\displaystyle \displaystyle {[-1,1]}}$ can be written ${\displaystyle \displaystyle {[\pi ,0]}}$

thus we have

${\displaystyle I=\int _{-1}^{+1}f(x)dx=\int _{\pi }^{0}-f(cos\theta )sin\theta d\theta =\int _{0}^{\pi }f(cos\theta )sin\theta d\theta }$

Author

--Egm6341.s10.team2.niki 18:47, 23 April 2010 (UTC)

Statement

Pg.42-3

We have the cosine series expressed as ${\displaystyle f(cos\theta )={\frac {a_{0}}{2}}+\sum _{k=1}^{\infty }a_{k}cos(k\theta )}$, we need to express teh constants ${\displaystyle a_{k}}$ as

${\displaystyle a_{k}={\frac {2}{\pi }}\int _{0}^{\pi }f(cos\theta )cos(k\theta )d\theta }$

Solution

We have the given expression for the cosine series as

${\displaystyle f(cos\theta )={\frac {a_{0}}{2}}+\sum _{k=1}^{\infty }a_{k}cos(k\theta )}$

multiplying both sides by ${\displaystyle cos(m\theta )}$ where k is not the same as m and integrating we get,

${\displaystyle \int _{0}^{2\pi }f(cos\theta )(cos(m\theta ))=\int _{0}^{2\pi }{\frac {a_{0}cos(m\theta )}{2}}+\int _{0}^{2\pi }\sum _{k=1}^{\infty }a_{k}cos(k\theta )cos(m\theta )}$

Using the property of orthogonality we know that ${\displaystyle \int _{0}^{2\pi }a_{k}cos(k\theta )cos(m\theta )}$ exists only when k = m i.e

${\displaystyle \int _{0}^{2\pi }f(cos\theta )(cos(k\theta ))d\theta ={\frac {a_{0}}{2}}\int _{0}^{2\pi }cos(m\theta )d\theta +a_{k}\int _{0}^{2\pi }cos^{2}(k\theta )d\theta }$

wkt,

${\displaystyle \int _{0}^{2\pi }cos(m\theta )d\theta =0}$ and

${\displaystyle \int _{0}^{2\pi }cos^{2}(k\theta )d\theta ={\frac {2\pi -0}{2}}=\pi }$

Thus we have by substituting and rearranging terms,

${\displaystyle a_{k}={\frac {2}{\pi }}\int _{0}^{\pi }f(cos\theta )cos(k\theta )d\theta }$

Author

--Egm6341.s10.team2.niki 14:43, 23 April 2010 (UTC)

Statement

Pg.42-3 Simplification of Clenshaw-Curtis Quadrature Rule

Solution

Author