Contents

Script file for Hw #8

% Here is the function of two variables whose  slope field we plot.
type fun1

function z = fun1(t,y)
z = t.^2-y;

Problem 1) Plotting the slope field and the four curves

close all
T=0:0.2:4;
Y=0:0.4:10;
[t,y]=meshgrid(T,Y);
 dt=ones(size(t)); dy=fun1(t,y);
 quiver(t,y,dt,dy)
 hold on
% Ploting
%
% $$y(t) =C e^{-t} + t^2 -2t + 2,  for C=-2, -1,  1, 2 $$
%
ym2=-2*exp(-T) + T.^2-2*T +2;
ym1=-1*exp(-T) + T.^2-2*T +2;
yp1=1*exp(-T) + T.^2-2*T +2;
yp2=2*exp(-T) + T.^2-2*T +2;
plot(T,ym2, T,ym1, T,yp1, T,yp2)

%

Problem 2) Find the exact solution of the IVP

.

%
% Imposing  y(0)=1  in y(t) =C e^{-t} + t^2 -2t + 2, we get C=-1.
% Thus the exact solution is y(t) =-e^{-t} + t^2 -2t + 2
%

Problem 3) Two steps of Euler's method

start with t0=0, w0=1, h=0.2. For i=0,1 DO

%  t_{i+1} = t_i +h,  w_{i+1} =w_i + h(t_i^2 -w_i)

% Computing  t, w using euler.m.
[t,w] = euler(@fun1,[0,0.4],1,0.2)
%

t =

         0    0.2000    0.4000


w =

    1.0000    0.8000    0.6480

Problem 4) Two steps of Heun's (Trapezoid) method

start with t0=0, w0=1, h=0.2. For i=0,1 DO

% t_{i+1} = t_i +h,
% f1 = f(t_i,w_i)   f2 = f(t_i+h, w_i +h*f1)$$
% w_{i+1} =w_i + h/2 (f1 + f2)

% Computing  t, w using heunT.m.
[t,w] = heunT(@fun1,[0,0.4],1,0.2)
%

t =

         0    0.2000    0.4000


w =

    1.0000    0.8240    0.6949

Problem 5) One steps of RK4 method

start with t0=0, w0=1, h=0.2. t_1 = t_0 +h,

% f1 = f(t_0,w_0); f2 = f(t_0+h/2, w_0 +h/2 *f1);
% f3 = f(t_0+h/2,w_0+h/2 *f2); f4 = f(t_0+h,w_i+h*f3);
%
%$$ w_1 =w_0 + h/6*(f1+2*f2+2*f3+f4) $$
% Computing  t, w using rk4.m.
[t,w] = rk4(@fun1,[0,0.2],1,0.2)
%

t =

         0    0.2000


w =

    1.0000    0.8213

Problem 6) The Resonant Spring model

Denote x'=y, Then x"=y' and we have THE IVP system

% x'=y
% y'=-25 x + 8 sin(5t)
% x(0) =0; y(0)=0.

% the function for ths RHS to be used as the first argument for rk4s.m is:
type VFrs

% To solve the system by rk4s, just do
[T, Z]=rk4s(@VFrs, [0, 2], [0,0],0.05);
Z'
%the first column of Z gives the numerical solution. To plot it DO
hold off
plot(T,Z(:,1))
%

function V=VFrs(t, z)
x=z(1);
y=z(2);
V=zeros(1,2);
V(1)=y;
V(2)=-25*x + 8*sin(5*t);
end

ans =

  Columns 1 through 14

         0    0.0008    0.0065    0.0213    0.0482    0.0888    0.1426    0.2073    0.2786    0.3506    0.4162    0.4677    0.4977    0.4996
         0    0.0495    0.1918    0.4090    0.6731    0.9489    1.1969    1.3775    1.4548    1.4005    1.1969    0.8397    0.3388   -0.2812

  Columns 15 through 28

    0.4683    0.4009    0.2972    0.1602   -0.0046   -0.1884   -0.3803   -0.5675   -0.7364   -0.8735   -0.9664   -1.0047   -0.9812   -0.8924
   -0.9820   -1.7145   -2.4215   -3.0427   -3.5189   -3.7971   -3.8355   -3.6075   -3.1044   -2.3383   -1.3415   -0.1663    1.1181    2.4296

  Columns 29 through 41

   -0.7393   -0.5271   -0.2660    0.0303    0.3444    0.6566    0.9463    1.1930    1.3778    1.4851    1.5036    1.4271    1.2555
    3.6785    4.7732    5.6274    6.1660    6.3316    6.0891    5.4298    4.3734    2.9679    1.2877   -0.5700   -2.4910   -4.3508

Problem 7) Rabits vs Foxes model

The Vector function implementing the RHS of the Rabit-Fox model is:

type VRF
% Solving  part a)
[T, Z]=rk4s(@VRF, [0, 5], [3000,120],0.1);
plot(Z(:,1),Z(:,2))
% Solving  part b)
[T, Z]=rk4s(@VRF, [0, 5], [5000,100],0.1);
hold on
plot(Z(:,1),Z(:,2))
% Comment: The model confirmes the ``CIRCLE'' (or the EGG) of LIFE!

function V=VRF(t, z)
x=z(1);
y=z(2);
V=zeros(1,2);
V(1)=2*x -0.02*x.*y;
V(2)=0.0002*x.*y -0.8*y;
end

Published with MATLAB® 7.8