Soil Temperature
Table of contents
Code
Main program
sp_05_02.m | View on GitHub
% Supplemental program 5.2
% ---------------------------------------------------------
% Diurnal cycle of soil temperature with phase change using
% "excess heat" or "apparent heat capacity"
% ---------------------------------------------------------
% --- Physical constants in physcon structure
physcon.tfrz = 273.15; % Freezing point of water (K)
physcon.cwat = 4188.0; % Specific heat of water (J/kg/K)
physcon.cice = 2117.27; % Specific heat of ice (J/kg/K)
physcon.rhowat = 1000.0; % Density of water (kg/m3)
physcon.rhoice = 917.0; % Density of ice (kg/m3)
physcon.cvwat = physcon.cwat * physcon.rhowat; % Heat capacity of water (J/m3/K)
physcon.cvice = physcon.cice * physcon.rhoice; % Heat capacity of ice (J/m3/K)
physcon.tkwat = 0.57; % Thermal conductivity of water (W/m/K)
physcon.tkice = 2.29; % Thermal conductivity of ice (W/m/K)
physcon.hfus = 0.3337e6; % Heat of fusion for water at 0 C (J/kg)
% --- Initialize soil texture variables
% Soil texture classes (Cosby et al. 1984. Water Resources Research 20:682-690)
% 1: sand
% 2: loamy sand
% 3: sandy loam
% 4: silty loam
% 5: loam
% 6: sandy clay loam
% 7 silty clay loam
% 8: clay loam
% 9: sandy clay
% 10: silty clay
% 11: clay
soilvar.silt = [ 5.0, 12.0, 32.0, 70.0, 39.0, 15.0, 56.0, 34.0, 6.0, 47.0, 20.0]; % Percent silt
soilvar.sand = [92.0, 82.0, 58.0, 17.0, 43.0, 58.0, 10.0, 32.0, 52.0, 6.0, 22.0]; % Percent sand
soilvar.clay = [ 3.0, 6.0, 10.0, 13.0, 18.0, 27.0, 34.0, 34.0, 42.0, 47.0, 58.0]; % Percent clay
% Volumetric soil water content at saturation (porosity)
% (Clapp and Hornberger. 1978. Water Resources Research 14:601-604)
soilvar.watsat = [0.395, 0.410, 0.435, 0.485, 0.451, 0.420, 0.477, 0.476, 0.426, 0.492, 0.482];
% --- Model run control parameters
tmean = physcon.tfrz + 15.0; % Mean daily air temperature for diurnal cycle (K)
trange = 10.0; % Temperature range for diurnal cycle (K)
dt = 1800; % Time step (seconds)
nday = 200; % Number of days
soilvar.soil_texture = 1; % Soil texture class: sand
soilvar.method = 'excess-heat'; % Use excess heat for phase change
%soilvar.soil_texture = 11; % Soil texture class: clay
%soilvar.method = 'apparent-heat-capacity'; % Use apparent heat capacity for phase change
% --- Initialize soil layer variables
% Number of layers in soil profile
soilvar.nsoi = 120;
% Soil layer thickness (m)
for i = 1:soilvar.nsoi
soilvar.dz(i) = 0.025;
end
% Soil depth (m) at i+1/2 interface between layers i and i+1 (negative distance from surface)
soilvar.z_plus_onehalf(1) = -soilvar.dz(1);
for i = 2:soilvar.nsoi
soilvar.z_plus_onehalf(i) = soilvar.z_plus_onehalf(i-1) - soilvar.dz(i);
end
% Soil depth (m) at center of layer i (negative distance from surface)
soilvar.z(1) = 0.5 * soilvar.z_plus_onehalf(1);
for i = 2:soilvar.nsoi
soilvar.z(i) = 0.5 * (soilvar.z_plus_onehalf(i-1) + soilvar.z_plus_onehalf(i));
end
% Thickness between between z(i) and z(i+1)
for i = 1:soilvar.nsoi-1
soilvar.dz_plus_onehalf(i) = soilvar.z(i) - soilvar.z(i+1);
end
soilvar.dz_plus_onehalf(soilvar.nsoi) = 0.5 * soilvar.dz(soilvar.nsoi);
% Initial soil temperature (K) and unfrozen and frozen water (kg H2O/m2)
for i = 1:soilvar.nsoi
% Temperature (K)
soilvar.tsoi(i) = physcon.tfrz + 2.0;
% Soil water at saturation (kg H2O/m2)
h2osoi_sat = soilvar.watsat(soilvar.soil_texture) * physcon.rhowat * soilvar.dz(i);
% Actual water content is some fraction of saturation
if (soilvar.tsoi(i) > physcon.tfrz)
soilvar.h2osoi_ice(i) = 0;
soilvar.h2osoi_liq(i) = 0.8 * h2osoi_sat;
else
soilvar.h2osoi_liq(i) = 0;
soilvar.h2osoi_ice(i) = 0.8 * h2osoi_sat;
end
end
% --- Time stepping loop to increment soil temperature
% Counter for output file
m = 0;
% Main loop is NTIM iterations per day with a time step of DT seconds.
% This is repeated NDAY times.
ntim = round(86400/dt);
for iday = 1:nday
fprintf('day = %6.0f\n',iday)
for itim = 1:ntim
% Hour of day
hour = itim * (dt/86400 * 24);
% Surface temperature: Constant value TMEAN if TRANGE = 0. Otherwise, use a sine
% wave with max (TMEAN + 1/2 TRANGE) at 2 pm and min (TMEAN - 1/2 TRANGE) at 2 am
tsurf = tmean + 0.5 * trange * sin(2*pi/24 * (hour-8.0));
% Thermal conductivity and heat capacity
[soilvar] = soil_thermal_properties (physcon, soilvar);
% Soil temperatures
[soilvar] = soil_temperature (physcon, soilvar, tsurf, dt);
% Save hourly soil temperature profile for last day
if (iday == nday)
% Surface output
% Vector format - to write to data file
m = m + 1;
hour_vec(m) = hour;
z_vec(m) = 0;
tsoi_vec(m) = tsurf - physcon.tfrz; % deg C
% For MATLAB contour
hour_out(itim) = hour;
z_out(1) = 0;
tsoi_out(1,itim) = tsurf - physcon.tfrz; % deg C
% Soil layers for top 100 cm
for i = 1:soilvar.nsoi
if (soilvar.z(i) > -1.0)
m = m + 1;
hour_vec(m) = hour;
z_vec(m) = soilvar.z(i) * 100; % cm
tsoi_vec(m) = soilvar.tsoi(i) - physcon.tfrz; % deg C
z_out(i+1) = soilvar.z(i) * 100; % cm
tsoi_out(i+1,itim) = soilvar.tsoi(i) - physcon.tfrz; % deg C
end
end
end
end
end
% --- Write output file
A = [hour_vec; z_vec; tsoi_vec];
fileID = fopen('data.txt','w');
fprintf(fileID,'%12s %12s %12s\n','hour','z','tsoi');
fprintf(fileID,'%12.3f %12.3f %12.3f\n', A);
fclose(fileID);
% --- Make contour plot
contour(hour_out,z_out,tsoi_out,'ShowText','on')
title('Soil Temperature (^oC)')
xlabel('Time of day (hours)')
ylabel('Soil depth (cm)')
Aux. programs
phase_change.m | View on GitHub
function [soilvar] = phase_change (physcon, soilvar, dt)
% Adjust temperatures for phase change. Freeze or melt ice using
% energy excess or deficit needed to change temperature to the
% freezing point.
% ------------------------------------------------------
% Input
% dt ! Time step (s)
% physcon.hfus ! Heat of fusion for water at 0 C (J/kg)
% physcon.tfrz ! Freezing point of water (K)
% soilvar.nsoi ! Number of soil layers
% soilvar.dz ! Soil layer thickness (m)
% soilvar.cv ! Volumetric heat capacity (J/m3/K)
%
% Input/output
% soilvar.tsoi ! Soil temperature (K)
% soilvar.h2osoi_liq ! Unfrozen water, liquid (kg H2O/m2)
% soilvar.h2osoi_ice ! Frozen water, ice (kg H2O/m2)
%
% Output
% soilvar.hfsoi ! Soil phase change energy flux (W/m2)
% ------------------------------------------------------
% --- Initialize total soil heat of fusion to zero
soilvar.hfsoi = 0;
% --- Now loop over all soil layers to calculate phase change
for i = 1:soilvar.nsoi
% --- Save variables prior to phase change
wliq0 = soilvar.h2osoi_liq(i); % Amount of liquid water before phase change
wice0 = soilvar.h2osoi_ice(i); % Amount of ice before phase change
wmass0 = wliq0 + wice0; % Amount of total water before phase change
tsoi0 = soilvar.tsoi(i); % Soil temperature before phase change
% --- Identify melting or freezing layers and set temperature to freezing
% Default condition is no phase change (imelt = 0)
imelt = 0;
% Melting: if ice exists above melt point, melt some to liquid.
% Identify melting by imelt = 1
if (soilvar.h2osoi_ice(i) > 0 & soilvar.tsoi(i) > physcon.tfrz)
imelt = 1;
soilvar.tsoi(i) = physcon.tfrz;
end
% Freezing: if liquid exists below melt point, freeze some to ice.
% Identify freezing by imelt = 2
if (soilvar.h2osoi_liq(i) > 0 & soilvar.tsoi(i) < physcon.tfrz)
imelt = 2;
soilvar.tsoi(i) = physcon.tfrz;
end
% --- Calculate energy for freezing or melting
% The energy for freezing or melting (W/m2) is assessed from the energy
% excess or deficit needed to change temperature to the freezing point.
% This is a potential energy flux, because cannot melt more ice than is
% present or freeze more liquid water than is present.
%
% heat_flux_pot > 0: freezing; heat_flux_pot < 0: melting
if (imelt > 0)
heat_flux_pot = (soilvar.tsoi(i) - tsoi0) * soilvar.cv(i) * soilvar.dz(i) / dt;
else
heat_flux_pot = 0;
end
% Maximum energy for melting or freezing (W/m2)
if (imelt == 1)
heat_flux_max = -soilvar.h2osoi_ice(i) * physcon.hfus / dt;
end
if (imelt == 2)
heat_flux_max = soilvar.h2osoi_liq(i) * physcon.hfus / dt;
end
% --- Now freeze or melt ice
if (imelt > 0)
% Change in ice (kg H2O/m2/s): freeze (+) or melt (-)
ice_flux = heat_flux_pot / physcon.hfus;
% Update ice (kg H2O/m2)
soilvar.h2osoi_ice(i) = wice0 + ice_flux * dt;
% Cannot melt more ice than is present
soilvar.h2osoi_ice(i) = max(0, soilvar.h2osoi_ice(i));
% Ice cannot exceed total water that is present
soilvar.h2osoi_ice(i) = min(wmass0, soilvar.h2osoi_ice(i));
% Update liquid water (kg H2O/m2) for change in ice
soilvar.h2osoi_liq(i) = max(0, (wmass0-soilvar.h2osoi_ice(i)));
% Actual energy flux from phase change (W/m2). This is equal to
% heat_flux_pot except if tried to melt too much ice.
heat_flux = physcon.hfus * (soilvar.h2osoi_ice(i) - wice0) / dt;
% Sum energy flux from phase change (W/m2)
soilvar.hfsoi = soilvar.hfsoi + heat_flux;
% Residual energy not used in phase change is added to soil temperature
residual = heat_flux_pot - heat_flux;
soilvar.tsoi(i) = soilvar.tsoi(i) - residual * dt / (soilvar.cv(i) * soilvar.dz(i));
% Error check: make sure actual phase change does not exceed permissible phase change
if (abs(heat_flux) > abs(heat_flux_max))
error ('Soil temperature energy conservation error: phase change')
end
% Freezing: make sure actual phase change does not exceed permissible phase change
% and that the change in ice does not exceed permissible change
if (imelt == 2)
% Energy flux (W/m2)
constraint = min(heat_flux_pot, heat_flux_max);
err = heat_flux - constraint;
if (abs(err) > 1e-03)
error ('Soil temperature energy conservation error: freezing energy flux')
end
% Change in ice (kg H2O/m2)
err = (soilvar.h2osoi_ice(i) - wice0) - constraint / physcon.hfus * dt;
if (abs(err) > 1e-03)
error ('Soil temperature energy conservation error: freezing ice flux')
end
end
% Thawing: make sure actual phase change does not exceed permissible phase change
% and that the change in ice does not exceed permissible change
if (imelt == 1)
% Energy flux (W/m2)
constraint = max(heat_flux_pot, heat_flux_max);
err = heat_flux - constraint;
if (abs(err) > 1e-03)
error ('Soil temperature energy conservation error: thawing energy flux')
end
% Change in ice (kg H2O/m2)
err = (soilvar.h2osoi_ice(i) - wice0) - constraint / physcon.hfus * dt;
if (abs(err) > 1e-03)
error ('Soil temperature energy conservation error: thawing ice flux')
end
end
end
end
soil_temperature.m | View on GitHub
function [soilvar] = soil_temperature (physcon, soilvar, tsurf, dt)
% Use an implicit formulation with the surface boundary condition specified
% as the surface temperature to solve for soil temperatures at time n+1.
%
% Calculate soil temperatures as:
%
% dT d dT
% cv -- = -- (k --)
% dt dz dz
%
% where: T = temperature (K)
% t = time (s)
% z = depth (m)
% cv = volumetric heat capacity (J/m3/K)
% k = thermal conductivity (W/m/K)
%
% Set up a tridiagonal system of equations to solve for T at time n+1,
% where the temperature equation for layer i is
%
% d_i = a_i [T_i-1] n+1 + b_i [T_i] n+1 + c_i [T_i+1] n+1
%
% For soil layers undergoing phase change, set T_i = Tf (freezing) and use
% excess energy to freeze or melt ice:
%
% Hf_i = (Tf - [T_i] n+1) * cv_i * dz_i / dt
%
% During the phase change, the unfrozen and frozen soil water
% (h2osoi_liq, h2osoi_ice) are adjusted.
%
% Or alternatively, use the apparent heat capacity method to
% account for phase change. In this approach, h2osoi_liq
% and h2osoi_ice are not calculated.
%
% ------------------------------------------------------
% Input
% tsurf ! Surface temperature (K)
% dt ! Time step (s)
% soilvar.method ! Use excess heat or apparent heat capacity for phase change
% soilvar.nsoi ! Number of soil layers
% soilvar.z ! Soil depth (m)
% soilvar.z_plus_onehalf ! Soil depth (m) at i+1/2 interface between layers i and i+1
% soilvar.dz ! Soil layer thickness (m)
% soilvar.dz_plus_onehalf ! Thickness (m) between between i and i+1
% soilvar.tk ! Thermal conductivity (W/m/K)
% soilvar.cv ! Heat capacity (J/m3/K)
%
% Input/output
% soilvar.tsoi ! Soil temperature (K)
% soilvar.h2osoi_liq ! Unfrozen water, liquid (kg H2O/m2)
% soilvar.h2osoi_ice ! Frozen water, ice (kg H2O/m2)
%
% Output
% soilvar.gsoi ! Energy flux into soil (W/m2)
% soilvar.hfsoi ! Soil phase change energy flux (W/m2)
% ------------------------------------------------------
% --- Save current soil temperature for energy conservation check
for i = 1:soilvar.nsoi
tsoi0(i) = soilvar.tsoi(i);
end
% --- Thermal conductivity at interface (W/m/K)
for i = 1:soilvar.nsoi-1
tk_plus_onehalf(i) = soilvar.tk(i) * soilvar.tk(i+1) * (soilvar.z(i)-soilvar.z(i+1)) / ...
(soilvar.tk(i)*(soilvar.z_plus_onehalf(i)-soilvar.z(i+1)) + soilvar.tk(i+1)*(soilvar.z(i)-soilvar.z_plus_onehalf(i)));
end
% --- Set up tridiagonal matrix
% Top soil layer with tsurf as boundary condition
i = 1;
m = soilvar.cv(i) * soilvar.dz(i) / dt;
a(i) = 0;
c(i) = -tk_plus_onehalf(i) / soilvar.dz_plus_onehalf(i);
b(i) = m - c(i) + soilvar.tk(i) / (0 - soilvar.z(i));
d(i) = m * soilvar.tsoi(i) + soilvar.tk(i) / (0 - soilvar.z(i)) * tsurf;
% Layers 2 to nsoi-1
for i = 2:soilvar.nsoi-1
m = soilvar.cv(i) * soilvar.dz(i) / dt;
a(i) = -tk_plus_onehalf(i-1) / soilvar.dz_plus_onehalf(i-1);
c(i) = -tk_plus_onehalf(i) / soilvar.dz_plus_onehalf(i);
b(i) = m - a(i) - c(i);
d(i) = m * soilvar.tsoi(i);
end
% Bottom soil layer with zero heat flux
i = soilvar.nsoi;
m = soilvar.cv(i) * soilvar.dz(i) / dt;
a(i) = -tk_plus_onehalf(i-1) / soilvar.dz_plus_onehalf(i-1);
c(i) = 0;
b(i) = m - a(i);
d(i) = m * soilvar.tsoi(i);
% --- Solve for soil temperature
[soilvar.tsoi] = tridiagonal_solver (a, b, c, d, soilvar.nsoi);
% --- Derive energy flux into soil (W/m2)
soilvar.gsoi = soilvar.tk(1) * (tsurf - soilvar.tsoi(1)) / (0 - soilvar.z(1));
% --- Phase change for soil layers undergoing freezing of thawing
switch soilvar.method
case 'apparent-heat-capacity'
% No explicit phase change energy flux. This is included in the heat capacity.
soilvar.hfsoi = 0;
case 'excess-heat'
% Adjust temperatures for phase change. Freeze or melt ice using energy
% excess or deficit needed to change temperature to the freezing point.
% The variable hfsoi is returned as the energy flux from phase change (W/m2).
[soilvar] = phase_change (physcon, soilvar, dt);
end
% --- Check for energy conservation
% Sum change in energy (W/m2)
edif = 0;
for i = 1:soilvar.nsoi
edif = edif + soilvar.cv(i) * soilvar.dz(i) * (soilvar.tsoi(i) - tsoi0(i)) / dt;
end
% Error check
err = edif - soilvar.gsoi - soilvar.hfsoi;
if (abs(err) > 1e-03)
error ('Soil temperature energy conservation error')
end
soil_thermal_properties.m | View on GitHub
function [soilvar] = soil_thermal_properties (physcon, soilvar)
% Calculate soil thermal conductivity and heat capacity
% ------------------------------------------------------
% Input
% physcon.hfus ! Heat of fusion for water at 0 C (J/kg)
% physcon.tfrz ! Freezing point of water (K)
% physcon.tkwat ! Thermal conductivity of water (W/m/K)
% physcon.tkice ! Thermal conductivity of ice (W/m/K)
% physcon.cvwat ! Heat capacity of water (J/m3/K)
% physcon.cvice ! Heat capacity of ice (J/m3/K)
% physcon.rhowat ! Density of water (kg/m3)
% physcon.rhoice ! Density of ice (kg/m3)
% soilvar.method ! Use excess heat or apparent heat capacity for phase change
% soilvar.soil_texture ! Soil texture class
% soilvar.sand ! Percent sand
% soilvar.watsat ! Volumetric soil water content at saturation (porosity)
% soilvar.nsoi ! Number of soil layers
% soilvar.dz ! Soil layer thickness (m)
% soilvar.tsoi ! Soil temperature (K)
% soilvar.h2osoi_liq ! Unfrozen water, liquid (kg H2O/m2)
% soilvar.h2osoi_ice ! Frozen water, ice (kg H2O/m2)
%
% Input/output
% soilvar.tk ! Thermal conducitivty (W/m/K)
% soilvar.cv ! Volumetric heat capacity (J/m3/K)
% ------------------------------------------------------
for i = 1:soilvar.nsoi
% --- Soil texture to process
k = soilvar.soil_texture;
% --- Volumetric soil water and ice
watliq = soilvar.h2osoi_liq(i) / (physcon.rhowat * soilvar.dz(i));
watice = soilvar.h2osoi_ice(i) / (physcon.rhoice * soilvar.dz(i));
% Fraction of total volume that is liquid water
fliq = watliq / (watliq + watice);
% Soil water relative to saturation
s = min((watliq + watice)/soilvar.watsat(k), 1);
% --- Dry thermal conductivity (W/m/K) from bulk density (kg/m3)
bd = 2700 * (1 - soilvar.watsat(k));
tkdry = (0.135 * bd + 64.7) / (2700 - 0.947 * bd);
% --- Soil solids thermal conducitivty (W/m/K)
% Thermal conductivity of quartz (W/m/K)
tk_quartz = 7.7;
% Quartz fraction
quartz = soilvar.sand(k) / 100;
% Thermal conductivity of other minerals (W/m/K)
if (quartz > 0.2)
tko = 2;
else
tko = 3;
end
% Thermal conductivity of soil solids (W/m/K)
tksol = tk_quartz^quartz * tko^(1-quartz);
% --- Saturated thermal conductivity (W/m/K) and unfrozen and frozen values
tksat = tksol^(1-soilvar.watsat(k)) * physcon.tkwat^(fliq*soilvar.watsat(k)) * ...
physcon.tkice^(soilvar.watsat(k)-fliq*soilvar.watsat(k));
tksat_u = tksol^(1-soilvar.watsat(k)) * physcon.tkwat^soilvar.watsat(k);
tksat_f = tksol^(1-soilvar.watsat(k)) * physcon.tkice^soilvar.watsat(k);
% --- Kersten number and unfrozen and frozen values
if (soilvar.sand(k) < 50)
ke_u = log10(max(s,0.1)) + 1;
else
ke_u = 0.7 * log10(max(s,0.05)) + 1;
end
ke_f = s;
if (soilvar.tsoi(i) >= physcon.tfrz)
ke = ke_u;
else
ke = ke_f;
end
% --- Thermal conductivity (W/m/K) and unfrozen and frozen values
soilvar.tk(i) = (tksat - tkdry) * ke + tkdry;
tku = (tksat_u - tkdry) * ke_u + tkdry;
tkf = (tksat_f - tkdry) * ke_f + tkdry;
% --- Heat capacity of soil solids (J/m3/K)
cvsol = 1.926e06;
% --- Heat capacity (J/m3/K) and unfrozen and frozen values
soilvar.cv(i) = (1 - soilvar.watsat(k)) * cvsol + physcon.cvwat * watliq + physcon.cvice * watice;
cvu = (1 - soilvar.watsat(k)) * cvsol + physcon.cvwat * (watliq + watice);
cvf = (1 - soilvar.watsat(k)) * cvsol + physcon.cvice * (watliq + watice);
% --- Adjust heat capacity and thermal conductivity if using apparent heat capacity
switch soilvar.method
case 'apparent-heat-capacity'
% Temperature range for freezing and thawing (K)
tinc = 0.5;
% Heat of fusion (J/m3) - This is equivalent to ql = hfus * (h2osoi_liq + h2osoi_ice) / dz
ql = physcon.hfus * (physcon.rhowat * watliq + physcon.rhoice * watice);
% Heat capacity and thermal conductivity
if (soilvar.tsoi(i) > physcon.tfrz+tinc)
soilvar.cv(i) = cvu;
soilvar.tk(i) = tku;
end
if (soilvar.tsoi(i) >= physcon.tfrz-tinc & soilvar.tsoi(i) <= physcon.tfrz+tinc)
soilvar.cv(i) = (cvf + cvu) / 2 + ql / (2 * tinc);
soilvar.tk(i) = tkf + (tku - tkf) * (soilvar.tsoi(i) - physcon.tfrz + tinc) / (2 * tinc);
end
if (soilvar.tsoi(i) < physcon.tfrz-tinc)
soilvar.cv(i) = cvf;
soilvar.tk(i) = tkf;
end
end
end
tridiagonal_solver.m | View on GitHub
function [u] = tridiagonal_solver (a, b, c, d, n)
% Solve for U given the set of equations R * U = D, where U is a vector
% of length N, D is a vector of length N, and R is an N x N tridiagonal
% matrix defined by the vectors A, B, C each of length N. A(1) and
% C(N) are undefined and are not referenced.
%
% |B(1) C(1) ... ... ... |
% |A(2) B(2) C(2) ... ... |
% R = | A(3) B(3) C(3) ... |
% | ... A(N-1) B(N-1) C(N-1)|
% | ... ... A(N) B(N) |
%
% The system of equations is written as:
%
% A_i * U_i-1 + B_i * U_i + C_i * U_i+1 = D_i
%
% for i = 1 to N. The solution is found by rewriting the
% equations so that:
%
% U_i = F_i - E_i * U_i+1
% --- Forward sweep (1 -> N) to get E and F
e(1) = c(1) / b(1);
for i = 2: 1: n-1
e(i) = c(i) / (b(i) - a(i) * e(i-1));
end
f(1) = d(1) / b(1);
for i = 2: 1: n
f(i) = (d(i) - a(i) * f(i-1)) / (b(i) - a(i) * e(i-1));
end
% --- Backward substitution (N -> 1) to solve for U
u(n) = f(n);
for i = n-1: -1: 1
u(i) = f(i) - e(i) * u(i+1);
end
Output
Figures
Figure 1
Text
data.txt | View on GitHub | View raw
sp_05_02_out.txt (standard output) | View on GitHub | View raw
day = 1
day = 2
day = 3
day = 4
day = 5
day = 6
day = 7
day = 8
day = 9
day = 10
day = 11
day = 12
day = 13
day = 14
day = 15
day = 16
day = 17
day = 18
day = 19
day = 20
day = 21
day = 22
day = 23
day = 24
day = 25
day = 26
day = 27
day = 28
day = 29
day = 30
day = 31
day = 32
day = 33
day = 34
day = 35
day = 36
day = 37
day = 38
day = 39
day = 40
day = 41
day = 42
day = 43
day = 44
day = 45
day = 46
day = 47
day = 48
day = 49
day = 50
day = 51
day = 52
day = 53
day = 54
day = 55
day = 56
day = 57
day = 58
day = 59
day = 60
day = 61
day = 62
day = 63
day = 64
day = 65
day = 66
day = 67
day = 68
day = 69
day = 70
day = 71
day = 72
day = 73
day = 74
day = 75
day = 76
day = 77
day = 78
day = 79
day = 80
day = 81
day = 82
day = 83
day = 84
day = 85
day = 86
day = 87
day = 88
day = 89
day = 90
day = 91
day = 92
day = 93
day = 94
day = 95
day = 96
day = 97
day = 98
day = 99
day = 100
day = 101
day = 102
day = 103
day = 104
day = 105
day = 106
day = 107
day = 108
day = 109
day = 110
day = 111
day = 112
day = 113
day = 114
day = 115
day = 116
day = 117
day = 118
day = 119
day = 120
day = 121
day = 122
day = 123
day = 124
day = 125
day = 126
day = 127
day = 128
day = 129
day = 130
day = 131
day = 132
day = 133
day = 134
day = 135
day = 136
day = 137
day = 138
day = 139
day = 140
day = 141
day = 142
day = 143
day = 144
day = 145
day = 146
day = 147
day = 148
day = 149
day = 150
day = 151
day = 152
day = 153
day = 154
day = 155
day = 156
day = 157
day = 158
day = 159
day = 160
day = 161
day = 162
day = 163
day = 164
day = 165
day = 166
day = 167
day = 168
day = 169
day = 170
day = 171
day = 172
day = 173
day = 174
day = 175
day = 176
day = 177
day = 178
day = 179
day = 180
day = 181
day = 182
day = 183
day = 184
day = 185
day = 186
day = 187
day = 188
day = 189
day = 190
day = 191
day = 192
day = 193
day = 194
day = 195
day = 196
day = 197
day = 198
day = 199
day = 200