Contents

KKKA05 Tenta 2010-03-11 CARMEN AREVALO

clear all
close all

Question 1

n=6;
for i=1:n
    for j=1:n
        A(i,j)=5/(i+2*j-1);
    end
end
A
x=ones(n,1);
b=A*x
xsol=A\b
res=b-A*xsol
infnorm=norm(res,inf)
A =
          2.5         1.25      0.83333        0.625          0.5      0.41667
       1.6667            1      0.71429      0.55556      0.45455      0.38462
         1.25      0.83333        0.625          0.5      0.41667      0.35714
            1      0.71429      0.55556      0.45455      0.38462      0.33333
      0.83333        0.625          0.5      0.41667      0.35714       0.3125
      0.71429      0.55556      0.45455      0.38462      0.33333      0.29412
b =
        6.125
       4.7757
       3.9821
       3.4423
       3.0446
       2.7365
xsol =
            1
            1
            1
            1
            1
            1
res =
 -8.8818e-016
            0
            0
            0
 -4.4409e-016
 -4.4409e-016
infnorm =
  8.8818e-016

Question 2

X=0.95;
g=@(r)log((1+r*(1-X))./(r*(1-X)))-(r+1)./(r.*(1+r*(1-X)));
r=linspace(0.1,1);
y=g(r);

A plot is drawn to get a reasonable guess to use with fzero

plot(r,y,r,0)

The optimal recycle ratio is

R=fzero(g,0.3)
R =
      0.30715

Question 3

p0=2555;
k=0.026;
pmax=12000;
pprime=@(t,p)k*(1-p/pmax).*p;
tdata=1950:10:2010;
[t,p] = ode45(pprime,tdata,p0);
format short g
sol=[t,p]
tspan=1950:10:2100;
options=odeset('events',@myevents10b);
[t,p,te,pe,ie] = ode45(pprime,tspan,p0,options);
ev=[te,pe]
sol =
         1950         2555
         1960       3116.6
         1970       3752.6
         1980       4453.4
         1990       5202.4
         2000       5977.7
         2010       6753.7
ev =
       2026.9         8000

The population will reach 8000 million in 2026

pdata=[2555 3040 3708 4454 5276 6079 6806];
figure(1)
plot(tdata,pdata,'o',t,p)
legend('data','simulation')
xlabel('year')
ylabel('million inhabitants')
title('Actual data and simulation of world population')

Event file

type myevents10b.m
function [value,isterminal,direction] = myevents10b(t,y)
%MYEVENTS10B Event function for prob 3
%   in tenta KKKA05 2010-03-11
value=y-8000;
isterminal=1;
direction=0;

end