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

(5)Use the code of HS algorithm and run the code for $h=2^{k}h$ and also develop and run the code for Forward and backward Eular methods for the same step size

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

For k = 2

For k=3

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

For k=2

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

Ref Lecture 42.2

Problem Statement

Prove that $a\int \limits _{t=0}^{2\pi }[1-e^{2}cos^{2}t]^{1/2}dt=4a\int \limits _{\alpha =0}^{\pi /2}[1-e^{2}sin^{2}\alpha ]^{1/2}d\alpha$

Solution

$a\int \limits _{t=0}^{2\pi }[1-e^{2}cos^{2}t]^{1/2}dt$

$=4a\int \limits _{t=0}^{\pi /2}[1-e^{2}cos^{2}t]^{1/2}dt$

Let $cost=Sin\alpha$

$-sintdt=cos\alpha d\alpha$

Using above values we obtain

$=4a(-)\int \limits _{\alpha =0}^{\pi /2}[1-e^{2}sin^{2}\alpha ]^{1/2}(-){\frac {cos\alpha }{sint}}d\alpha$

$=4a(-)\int \limits _{\alpha =0}^{\pi /2}[1-e^{2}sin^{2}\alpha ]^{1/2}(-){\frac {cos\alpha }{sint}}d\alpha$

$=4a\int \limits _{\alpha =0}^{\pi /2}[1-e^{2}sin^{2}\alpha ]^{1/2}{\frac {cos\alpha }{(1-sin^{2}\alpha )^{(}1/2)}}d\alpha$

$=4a\int \limits _{\alpha =0}^{\pi /2}[1-e^{2}sin^{2}\alpha ]^{1/2}{\frac {cos\alpha }{cos\alpha }}d\alpha$

$=4a\int \limits _{\alpha =0}^{\pi /2}[1-e^{2}sin^{2}\alpha ]^{1/2}d\alpha$

Hence Proved