University of Florida/Egm6341/s10.team3.aks/HW7

Ref Lecture p.40.2

Problem Statement

Integrate the logistic equation and use the code of question 3 for HS and run the code for $h=2^{k}h$ and also develop Forward Eular and Back Eular methods for the same step sizes.

Solution

Matlab code with change in step size is given below.

MATLAB CODE:

clc;
clear all;
h=0.01;      %Step-size
r=1.2;
xmax=10;
t=0:0.01:10; %Range of time,t 
x_0 = 2;     %Initial Condition
tol = 1e-6;  %Absolute Tolerance
syms x;
syms y;
f = r*x*(1-x/xmax);  %Given Verhulst function
 
for i=1:length(t)-1
    x(1)=x_0;
    x(i+1)=y;
    f1 = subs(f,x(i));      %The value of function at x(i) is calculated - Numerical value- Known
    f2 = subs(f,x(i+1));    %The value of function at x(i+1) is calculated - Numerical value- Unknown - A function of x(i+1)
 
    x_half = 0.5*(x(i)+x(i+1)) + (h* 2^(i-1)/8)*(f1-f2);   
    f_half = subs(f,x_half); %The value of function at x(i+ 1/2) is calculated - Unknown - A function of x(i+1)
 
    F = -x(i+1) + x(i) + (h* 2^(i-1)/6)*(f1 + f2 + 4*f_half); %The non-linear algebraic function F(x+1) = 0
 
    z_guess = x(i);  %Initial Guess
    flag = 0;
 
    while flag ~= 1
        new_z = z_guess - (1/(subs(diff(F),z_guess))) * subs(F,z_guess);  %Hermite-Simpson Rule applied
        if double(abs(new_z - z_guess)) <= tol                            %Check for absolute tolerance
             x(i+1) = double(new_z);  %If satisfied, this is the value of x(i+1) we are looking for
             flag = 1;
             break;
        else
            z_guess = new_z;        %If not satisfied, repeat iteration, by making the new guess of x(i+1) to be equal to the value obtained above.
            flag = 0;
        end
    end
end
 
plot(t,double(x));                  %Plot the function
xlabel('Time(t)')
ylabel('Population(x)');
title('Population vs Time');
grid on;

For k = 1

File:For k=1.png

For k = 2

File:For k=2.png

For k=3

File:For k=3.png

We can see that as we increase the value of k due to decrease in step size our curve shifted to the left side.

Below is the code for Forward Eular method

MATLAB CODE:

clc;
clear all;
h=0.01;      %Step-size
r=1.2;
xmax=10;
t=0:0.01:10; %Range of time,t 
x_0 = 2;     %Initial Condition
tol = 1e-6;  %Absolute Tolerance
syms x;
syms y;
f = r*x*(1-x/xmax);  %Given Verhulst function
 
for i=1:length(t)-1
    x(i)=x_0;
    x(i+1)=y;
    f1 = subs(f,x(i));      %The value of function at x(i) is calculated - Numerical value- Known
    f2 = subs(f,x(i+1));    %The value of function at x(i+1) is calculated - Numerical value- Unknown - A function of x(i+1)
 
    x_half = 0.5*(x(i)+x(i+1)) + (h/8)*(f1-f2);   
    f_half = subs(f,x_half); %The value of function at x(i+ 1/2) is calculated - Unknown - A function of x(i+1)
 
    F = -x(i+1) + x(i) + (h/6)*(f1 + f2 + 4*f_half); %The non-linear algebraic function F(x+1) = 0
 
    z_guess = x(i);  %Initial Guess
    flag = 0;
 
    while flag ~= 1
        new_z = z_guess * (1+ h*r);  % Forward Eular Method
        if double(abs(new_z - z_guess)) <= tol                            %Check for absolute tolerance
             x(i+1) = double(new_z);  %If satisfied, this is the value of x(i+1) we are looking for
             flag = 1;
             break;
        else
            z_guess = new_z;        %If not satisfied, repeat iteration, by making the new guess of x(i+1) to be equal to the value obtained above.
            flag = 0;
        end
    end
end
 
plot(t,double(x));                  %Plot the function
xlabel('Time(t)');
ylabel('Population(x)');
title('Population vs Time');
grid on;

Forward Eular breaks down at higher values of step size h .This shows that it is unstable in nature. Below is the code for Backword Eular Method

MATLAB CODE:

clc;
clear all;
h=0.01;      %Step-size
r=1.2;
xmax=10;
t=0:0.1:10; %Range of time,t 
x_0 = 2;     %Initial Condition
tol = 1e-6;  %Absolute Tolerance
syms x;
syms y;
f = r*x*(1-x/xmax);  %Given Verhulst function
 
for i=1:length(t)-1
    x(1)=x_0;
    x(i+1)=y;
    f1 = subs(f,x(i));      %The value of function at x(i) is calculated - Numerical value- Known
    f2 = subs(f,x(i+1));    %The value of function at x(i+1) is calculated - Numerical value- Unknown - A function of x(i+1)
 
    x_half = 0.5*(x(i)+x(i+1)) + (h/8)*(f1-f2);   
    f_half = subs(f,x_half); %The value of function at x(i+ 1/2) is calculated - Unknown - A function of x(i+1)
 
    F = -x(i+1) + x(i) + (h/6)*(f1 + f2 + 4*f_half); %The non-linear algebraic function F(x+1) = 0
 
    z_guess = x(i);  %Initial Guess
    flag = 0;
 
    while flag ~= 1
        new_z = z_guess / (1- h*r);  % backward Eular Method
        if double(abs(new_z - z_guess)) <= tol                            %Check for absolute tolerance
             x(i+1) = double(new_z);  %If satisfied, this is the value of x(i+1) we are looking for
             flag = 1;
             break;
        else
            z_guess = new_z;        %If not satisfied, repeat iteration, by making the new guess of x(i+1) to be equal to the value obtained above.
            flag = 0;
        end
    end
end
 
plot(t,double(x));                  %Plot the function
xlabel('Time(t)');
ylabel('Population(x)');
title('Population vs Time');
grid on;

For k = 1

File:BE k=1.png

For k=2

File:BE k=2.png

(12)Convert circumference of ellipse formula from t to α variable

Ref Lecture 42.2

Problem Statement

Prove that $a\underset{t=0}{\overset{2\pi }{\int }}[1-e^{2}co{s}^{2}t{]}^{1/2}dt=4a\underset{\alpha =0}{\overset{\pi /2}{\int }}[1-e^{2}si{n}^2\alpha {]}^{1/2}d\alpha$

Solution

$a\underset{t=0}{\overset{2\pi }{\int }}[1-e^{2}co{s}^{2}t{]}^{1/2}dt$

$=4a\underset{t=0}{\overset{\pi /2}{\int }}[1-e^{2}co{s}^{2}t{]}^{1/2}dt$

Let $cost=Sin\alpha$

$-sintdt=cos\alpha d\alpha$

Using above values we obtain

$=4a(-)\underset{\alpha =0}{\overset{\pi /2}{\int }}[1-e^{2}si{n}^2\alpha {]}^{1/2}(-)\frac{cos\alpha }{sint}d\alpha$

$=4a(-)\underset{\alpha =0}{\overset{\pi /2}{\int }}[1-e^{2}si{n}^2\alpha {]}^{1/2}(-)\frac{cos\alpha }{sint}d\alpha$

$=4a\underset{\alpha =0}{\overset{\pi /2}{\int }}[1-e^{2}si{n}^2\alpha {]}^{1/2}\frac{cos\alpha }{(1-si{n}^2\alpha {)}^{(}1/2)}d\alpha$

$=4a\underset{\alpha =0}{\overset{\pi /2}{\int }}[1-e^{2}si{n}^2\alpha {]}^{1/2}\frac{cos\alpha }{cos\alpha }d\alpha$

$=4a\underset{\alpha =0}{\overset{\pi /2}{\int }}[1-e^{2}si{n}^2\alpha {]}^{1/2}d\alpha$

Hence Proved