Jump to content

University of Florida/Egm6341/s10.Team2/HW7

From Wikiversity

problem 1: Stable Growth for Verhulst Equation

[edit | edit source]

Statement

[edit | edit source]

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

Solution

[edit | edit source]

Author

[edit | edit source]

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

problem 2: Solution of the Logistic equation

[edit | edit source]

Statement

[edit | edit source]

P.39-1

Solution

[edit | edit source]

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

     

Separating variables and integrating we have

     

which can be written as

     

We get

     
     

Rearranging we have

     

Author

[edit | edit source]

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


problem 3: Hermite-Simpson Algorithm to solve for the non-linear first order differential equation i.e. Verhulst population growth model on Pg.38-3

[edit | edit source]

Statement

[edit | edit source]

Pg.39-1

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

Consider the following: and
Case 1:
Case 2:

And also plot the two cases using matlab.

Solution

[edit | edit source]

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

[edit | edit source]

Solved by--Srikanth Madala (SM)


Problem4: Creating graphs, showing the bifurcation diagram for the logistic map and its sensitivity to initial conditions

[edit | edit source]

Statement

[edit | edit source]

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

[edit | edit source]
Creating Figure 15.6:Bifurcation Diagram
[edit | edit source]

Logistic Map:

The Bifurcation Diagram shows the period doubling process of the logistic map. As the constant r increases the period is doubled to 2 periods and then stability is lost which then leads to a doubling of period to 4 and so on.

The plot is as follows:

Matlab Code
[edit | edit source]

Matlab code used to create the graph as outlined in Pg. 455 of King et. al

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
Creating Figure 15.7: Sensitivity to Initial Conditions
[edit | edit source]

The Following Plot explores the big effect of a small change in the initial conditions. One can note the changes as 'n' becomes larger than 50. The graph was produced for the following Conditions:

The Following Equation was used to produce the iterations for the graph:

Matlab Code
[edit | edit source]

Matlab Code used to produce the aforementioned plot.

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

[edit | edit source]

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


problem 5: Variation of the plot with the increase in the time step size

[edit | edit source]

Statement

[edit | edit source]

Pg.40-2

Solution

[edit | edit source]

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

color scheme in the plot below
Color h-value
Red h=0.2
Blue h=0.4
Green h=0.8
Magenta h=1.6
Brown h=3.2

Author

[edit | edit source]

Solved by--Srikanth Madala (SM)


Problem6: Solving Ordinary Differential Equations using Euler Implicit Methods

[edit | edit source]

Statement

[edit | edit source]

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:


Initial Conditions are as follows:

The Physical Modeling Parameters to be used are












The Inputs are defined as follows:

Input T:

for
for
for

Input :

for
for as it varies linearly
for as it varies linearly
for as it varies linearly

Solution

[edit | edit source]
Defining the Inputs
[edit | edit source]

The first part of the solution is to create a function that will calculate the thrust and angle of a attack at a particular time. The following Matlab Codes were used:

Thrust:

% 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

Angle of Attack:

% 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
Technique used to solve the problem
[edit | edit source]

The overall goal is to find a solution for each of the states being considered (x,y,V,gamma).

From P 36.3, of the lectures, the following was established using the Hermite-Simpson Algorithm:

    

The Following is recognized and derived as:

    

It is now necessary to define the following using the previously mentioned Eq. 1:

    

The next step is to apply the Newton Raphson Method to solve the function

A Matlab Algorithm is used to apply the Newton Raphson Method as follows:

Newton Raphson Method
[edit | edit source]

Function for obtaining

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 for obtaining

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

Function for using the Newton Raphson Operations until an error threshold is reached

This Matlab code was developed using the published code of Team 3 as part of their HW 7 Report( HW 7: Team 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;
Results of Methods Used and comparison to ODE45
[edit | edit source]

Plot of Horizontal Position Vs. Time

Plot of Vertical Position(Height) Vs. Time

Plot of Velocity Vs. Time

Plot of Gamma Vs. Time

From each of these plots it can be concluded that the combination of the Hermite Simpson Algorith and the Newton Raphson Method is an effective way to solve the set of ordinary differential equations without much error, as compared to ODE45 solver of Matlab. The downside to this method is that it takes a high amount of computational resources and time.

Author

[edit | edit source]

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


problem7 :Integration of Logistic equation by

[edit | edit source]

Statement

[edit | edit source]

Integrate the Logistic equation using the inconsistent Trapezoidal-Simpson approximation

Pg.41-2

Solution

[edit | edit source]

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


for

Author

[edit | edit source]

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


problem8

[edit | edit source]

Statement

[edit | edit source]

Pg.42-1

Solution

[edit | edit source]

Author

[edit | edit source]

Problem9: Eliminating t from the Parameterization of an ellipse

[edit | edit source]

Statement

[edit | edit source]

Pg.42-1

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

Solution

[edit | edit source]






At this point t can be eliminated as a parameter.

Author

[edit | edit source]

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


problem10

[edit | edit source]

Statement

[edit | edit source]

Pg.42-2

Solution

[edit | edit source]

Author

[edit | edit source]

problem11: To show the integral for elliptical circumference

[edit | edit source]

Statement

[edit | edit source]

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

Solution

[edit | edit source]

Author

[edit | edit source]

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

problem12

[edit | edit source]

Statement

[edit | edit source]

Pg.42-2

Solution

[edit | edit source]

Author

[edit | edit source]

problem 13: Change of Variable

[edit | edit source]

Statement

[edit | edit source]

Pg.42-2

To show that

Solution

[edit | edit source]

We have the integral as

.

Using the transformation we have

and using the inverse transformation the limits can be written

thus we have

Author

[edit | edit source]

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


problem 14: Constants of the Cosine Series

[edit | edit source]

Statement

[edit | edit source]

Pg.42-3

We have the cosine series expressed as , we need to express teh constants as

Solution

[edit | edit source]

We have the given expression for the cosine series as

multiplying both sides by where k is not the same as m and integrating we get,

Using the property of orthogonality we know that exists only when k = m i.e

wkt,

and

Thus we have by substituting and rearranging terms,

Author

[edit | edit source]

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


problem15: Expression for Clenshaw-Curtis Quadrature Rule

[edit | edit source]

Statement

[edit | edit source]

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

Solution

[edit | edit source]

Author

[edit | edit source]

Signatures

[edit | edit source]

solutions of problems 2,7,13,14 --Egm6341.s10.team2.niki 18:48, 23 April 2010 (UTC)

solutions of problems 4,6, 9 --Egm6341.s10.Team2.GV 20:58, 23 April 2010 (UTC)

solutions of problems 1,11,15 --Egm6341.s10.team2.lee 22:47, 23 April 2010 (UTC)

Solution Authors and Reviewers

[edit | edit source]
HW Assignment # 7
Problem # Solution by Reviewed by
Problem 1 JP NN
Problem 2 NN SM
Problem 3 SM GV
Problem 4 GV JP
Problem 5 SM PO
Problem 6 GV JP
Problem 7 NN & SM GV
Problem 8 PO NN
Problem 9 GV PO
Problem 10 PO SM
Problem 11 JP PO
Problem 12 PO GV
Problem 13 NN NN
Problem 14 NN JP
Problem 15 JP SM