Mathematical Computing Assignment 2: Global Warming Model

This document presents a MATLAB solution to a global warming assignment. The assignment involves modeling the dynamics of carbon dioxide emissions and absorption in the Earth's atmosphere, shallow ocean, and deep ocean. The solution includes MATLAB code to solve a system of differential equations that model five key variables: the partial pressure of carbon dioxide in the atmosphere (p), dissolved carbon concentrations in the shallow (σs) and deep ocean (σd), and alkalinity in the shallow (αs) and deep ocean (αd). The code utilizes functions and solvers to analyze the rate of carbon emission and solve the differential equations, generating plots to visualize the carbon concentrations over time. The analysis includes an examination of how different parameters impact the CO2 levels in different regimes and the percentage of CO2 over time. The results show the changes in carbon concentrations from the year 1000 to 5000, with the maximum value of CO2 percentage in the year 5000. The solution aims to provide a simplified model for understanding the effect of burning fossil fuels on the Earth's atmosphere.

Global Warming
Name of the Student
Name of the University
Author Note

Question 1:
MATLAB code:
function ydot= fossil_fuels(t,y)
p=y(1); sigmas=y(2); sigmad=y(3); alphas=y(4); alphad=y(5);
r=[0.0 0.0 1.0 4.0 5.0 8.0 10.0 10.5 10.0 8.0 3.5 2.0 0.0 0.0];
yr=[1000 1850 1950 1980 2000 2050 2080 2100 2120 2150 2225 2300 2500 5000];
f=pchip(yr,r);
f = ppval(f,t);
d= 8.65;
miu1 = 4.95e2;
miu2 = 4.95e-2;
vs = 0.12;
vd = 1.23;
w = 1e-3;
k1 = 2.19e-4;
k2 = 6.12e-5;
k3 = 0.997148;
k4 = 6.79e-2;
hs = (sigmas - (sigmas^2 - k3*alphas*(2*sigmas - alphas))^(1/2))/k3;
cs = (alphas - hs)/2;

ps = k4*(hs^2)/cs;
ydot = [(ps - p)/d + f(end)/miu1,(1/vs)*((sigmad - sigmas)*w - k1 -
(ps-p)*miu2/d),(1/vd)*(k1 - (sigmad - sigmas)*w),(1/vs)*((alphad - alphas)*w -
k2),(1/vd)*(k2 - (alphad - alphas)*w)];
v = ydot;
u = linspace(1000,5000,length(ydot));
plot_handles=plot(u,v);
legend(plot_handles,'fossil fuel');
title('Rate of Emission of Carbon');
ylabel('rate');
xlabel('year');
end
Execution and Output:
>> t = [1000 1850 1950 1980 2000 2050 2080 2100 2120 2150 2225 2300 2500 5000];
>> y = [0.0 0.0 1.0 4.0 5.0 8.0 10.0 10.5 10.0 8.0 3.5 2.0 0.0 0.0];
>> ydot= fossil_fuels(t,y)
ydot =

Columns 1 through 11
0.0315 0.0315 0.0335 0.0395 0.0416 0.0476 0.0517 0.0527 0.0517 0.0476
0.0385
Columns 12 through 18
0.0355 0.0315 0.0315 -0.0065 -0.0006 0.0078 -0.0008
1000 1500 2000 2500 3000 3500 4000 4500 5000
year
-0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
rate
Rate of Emission of Carbon
fossil fuel
2. The given differential equations are

The differential equations are solved using the solvers available in MATLAB.
MATLAB code:
%% defining parameters
r=[0.0 0.0 1.0 4.0 5.0 8.0 10.0 10.5 10.0 8.0 3.5 2.0 0.0 0.0];
p=1.00; sigmas=2.01; sigmad=2.23; alphas=2.20; alphad=2.26;
yr=[1000 1850 1950 1980 2000 2050 2080 2100 2120 2150 2225 2300 2500 5000];
t = yr;
f=pchip(yr,r);
f = ppval(f,t);

d= 8.65;
miu1 = 4.95e2;
miu2 = 4.95e-2;
vs = 0.12;
vd = 1.23;
w = 1e-3;
k1 = 2.19e-4;
k2 = 6.12e-5;
k3 = 0.997148;
k4 = 6.79e-2;
hs = (sigmas - (sigmas^2 - k3*alphas*(2*sigmas - alphas))^(1/2))/k3;
cs = (alphas - hs)/2;
ps = k4*(hs^2)/cs;
%% Solving by ode45
intcon = [p,sigmas,sigmad,alphas,alphad]; % initial conditions
% variables assumed for defining odefunction : y(1) = p, y(2) = sigmas, y(3)= sigmad, y(4) =
alphas, y(5) = alphad

odefunc = @(t,y) [(ps - y(1))/d + f(4)/miu1; (1/vs)*((y(3) - y(2))*w - k1 - ((ps -
y(1))*miu2)/d);(1/vd)*(k1 - (y(3)-y(2))*w); (1/vs)*((y(5) - y(4))*w - k2); (1/vd)*(k2 - (y(5) -
y(4))*w)];
opts=odeset('reltol',1.e-6,'abstol',1.e-6,'stats','on');
[year,ycomp] = ode45(odefunc,[1000 5000],intcon,opts);
Output:
solvediffequs
154 successful steps
11 failed attempts
991 function evaluations
>> ycomp
ycomp =
1.0000 2.0100 2.2300 2.2000 2.2600
1.0121 2.0073 2.2300 2.2000 2.2600
1.0233 2.0050 2.2300 2.2000 2.2600
1.0337 2.0031 2.2300 2.2000 2.2600

1.0432 2.0015 2.2300 2.2000 2.2600
1.0515 2.0003 2.2300 2.2000 2.2600
1.0592 1.9993 2.2300 2.2000 2.2600
1.0663 1.9986 2.2300 2.2000 2.2600
1.0730 1.9981 2.2300 2.1999 2.2600
1.0795 1.9977 2.2300 2.1999 2.2600
1.0856 1.9976 2.2299 2.1999 2.2600
1.0912 1.9977 2.2299 2.1999 2.2600
1.0964 1.9979 2.2299 2.1999 2.2600
1.1015 1.9984 2.2299 2.1999 2.2600
1.1061 1.9990 2.2299 2.1999 2.2600
1.1104 1.9997 2.2299 2.1999 2.2600
1.1144 2.0006 2.2299 2.1999 2.2600
1.1182 2.0017 2.2299 2.1999 2.2600
1.1217 2.0029 2.2299 2.1999 2.2600
1.1249 2.0043 2.2299 2.1999 2.2600
1.1278 2.0057 2.2299 2.1999 2.2600
1.1307 2.0074 2.2299 2.1999 2.2600
1.1333 2.0092 2.2299 2.1999 2.2600
1.1356 2.0110 2.2299 2.1998 2.2600

1.1378 2.0129 2.2299 2.1998 2.2600
1.1398 2.0151 2.2299 2.1998 2.2600
1.1417 2.0173 2.2299 2.1998 2.2600
1.1434 2.0196 2.2299 2.1998 2.2600
1.1449 2.0219 2.2299 2.1998 2.2600
1.1463 2.0244 2.2299 2.1998 2.2600
1.1477 2.0270 2.2299 2.1998 2.2600
1.1488 2.0297 2.2299 2.1998 2.2600
1.1499 2.0323 2.2300 2.1998 2.2600
1.1509 2.0353 2.2300 2.1998 2.2600
1.1518 2.0382 2.2300 2.1998 2.2600
1.1526 2.0412 2.2300 2.1998 2.2600
1.1533 2.0441 2.2300 2.1997 2.2600
1.1539 2.0474 2.2301 2.1997 2.2600
1.1545 2.0506 2.2301 2.1997 2.2600
1.1550 2.0539 2.2302 2.1997 2.2600
1.1555 2.0572 2.2302 2.1997 2.2600
1.1559 2.0607 2.2302 2.1997 2.2600
1.1563 2.0643 2.2303 2.1997 2.2600
1.1566 2.0678 2.2304 2.1997 2.2600

1.1569 2.0713 2.2304 2.1997 2.2600
1.1571 2.0752 2.2305 2.1997 2.2600
1.1573 2.0791 2.2306 2.1997 2.2600
1.1575 2.0829 2.2306 2.1996 2.2600
1.1577 2.0867 2.2307 2.1996 2.2600
1.1578 2.0909 2.2308 2.1996 2.2600
The first few rows of the y components are displayed to suppress report length.
3.
Now, plot of each of the 5 y components are created with respect to values of time in years in
between 1000 and 5000.
MATLAB code for plot:
figure
plot_handle=plot(year,ycomp);
legend(plot_handle,'p','\sigma_s','\sigma_d','\alpha_s','\alpha_d','location','best');
title('Carbon in the atmosphere and ocean');
ylabel('carbon');
xlabel('year');
Plot:

1000 1500 2000 2500 3000 3500 4000 4500 5000
year
1
1.5
2
2.5
3
3.5
4
carbon
Carbon in the atmosphere and ocean
p
s
d
s
d
4. It can be seen from the above graph that the amounts of the carbon in regime p, αd and αs
are approximately same in year 5000 as they were in year 1000, whereas the amounts of
carbon differs by 1.5 ppm for the regimes σ s and σ d.
5. The percentage of C O2 reaches its maximum value in the year 5000 from the above graph.