Sharedgreenhouse.sagewsOpen in CoCalc
from scipy.constants import sigma     #importamos la constante de Stefan-Boltzmann


S = 250                    #Flujo solar W/m²
hf = 23.5                  #Coeficiente de convección forzada (velocidad del viento 5,36 m/s)
hc = 1.78                  #Coeficiente de convección libre turbulenta
kd = 200                   #Conductividad/distancia K/d
Ta = 288                   #Temperatura del medio
Tsky = 288                 #Temperatura característica de la radiación incidente desde arriba
epsilon= 1                 #emisividad absorción total
var('Ti,Te,TR')

f1(Ti, Te, TR) = sigma*epsilon*(Te^4-Tsky^4)+hf*(Te-Ta)-S
f2(Ti, Te, TR) = kd*(Ti-Te)-S
f3(Ti, Te, TR) = sigma*epsilon*(TR^4-Ti^4)+hc*(TR-Ti)^(4/3)-S

f(Ti, Te, TR) = (f1(Ti, Te, TR), f2(Ti, Te, TR), f3(Ti, Te, TR))

j = jacobian(f, [Ti, Te, TR])

x = vector([300, 300, 300]) # Initial values

data = [[0, x, n(norm(f(x[0], x[1], x[2])), digits=4)]]

for i in range(1,12):
    x = vector((n(d) for d in x - j(Ti=x[0], Te=x[1],
    TR=x[2]).inverse()*f(x[0], x[1], x[2])))
    data.append([i, x, norm(f(x[0], x[1], x[2]))])

# Table 

print "conveccion", hc*(data[11][1][2]-data[11][1][0])^(4/3),hf*(data[11][1][1]-Ta)
print "radiation", sigma*epsilon*((data[11][1][2])^4-(data[11][1][0])^4),sigma*epsilon*((data[11][1][1])^4-Tsky^4)

table([(data[i][0], data[i][1].n(digits=10),
           n(data[i][2], digits=4)) for i in range(0,12)],
           header_row = ["$i$", "$(T_i,T_e,T_R)$", "$norm(f)$"])
(Ti, Te, TR) conveccion 106.478122092393 201.441111513647 radiation 143.521877907607 48.5588884863539 $i$ $(T_i,T_e,T_R)$ $norm(f)$ +-----+-----------------------------------------+------------+ 0 (300.0000000, 300.0000000, 300.0000000) 367.8 1 (297.8340165, 296.5840165, 338.6570322) 299.9 2 (297.8219623, 296.5719623, 320.9880085) 23.40 3 (297.8219622, 296.5719622, 319.3476957) 0.2271 4 (297.8219622, 296.5719622, 319.3314623) 0.00002262 5 (297.8219622, 296.5719622, 319.3314607) 8.602e-13 6 (297.8219622, 296.5719622, 319.3314607) 1.256e-12 7 (297.8219622, 296.5719622, 319.3314607) 8.602e-13 8 (297.8219622, 296.5719622, 319.3314607) 1.256e-12 9 (297.8219622, 296.5719622, 319.3314607) 8.602e-13 10 (297.8219622, 296.5719622, 319.3314607) 1.256e-12 11 (297.8219622, 296.5719622, 319.3314607) 8.602e-13
from scipy.constants import sigma     #importamos la constante de Stefan-Boltzmann


S = 250                    #Flujo solar W/m²
hf = 23.5                  #Coeficiente de convección forzada (velocidad del viento 5,36 m/s)
hc = 1.78                  #Coeficiente de convección libre turbulenta
kd = 200                   #Conductividad/distancia K/d
Ta = 288                   #Temperatura del medio
Tsky = 288                 #Temperatura característica de la radiación incidente desde arriba
epsilon= 1                 #emisividad absorción total

var('Ti,Te,TR')

f1(Ti, Te, TR) = sigma*epsilon*(TR^4-Tsky^4)+hf*(Te-Ta)-S
f2(Ti, Te, TR) = sigma*epsilon*(TR^4-Tsky^4)+kd*(Ti-Te)-S
f3(Ti, Te, TR) = sigma*epsilon*(TR^4-Tsky^4)+hc*(TR-Ti)^(4/3)-S

f(Ti, Te, TR) = (f1(Ti, Te, TR), f2(Ti, Te, TR), f3(Ti, Te, TR))

j = jacobian(f, [Ti, Te, TR])

x = vector([300, 300, 300]) # Initial values

data = [[0, x, n(norm(f(x[0], x[1], x[2])), digits=4)]]

for i in range(1,12):
    x = vector((n(d) for d in x - j(Ti=x[0], Te=x[1],
    TR=x[2]).inverse()*f(x[0], x[1], x[2])))
    data.append([i, x, norm(f(x[0], x[1], x[2]))])

# Table 

print "conveccion", hc*(data[11][1][2]-data[11][1][0])^(4/3),hf*(data[11][1][1]-Ta)
print "radiation", sigma*epsilon*((data[11][1][2])^4-(Tsky)^4),sigma*epsilon*((data[11][1][2])^4-(Tsky)^4)

table([(data[i][0], data[i][1].n(digits=10),
           n(data[i][2], digits=4)) for i in range(0,12)],
           header_row = ["$i$", "$(T_i,T_e,T_R)$", "$norm(f)$"])
(Ti, Te, TR) conveccion 97.4830588020287 97.4830588020289 radiation 152.516941197971 152.516941197971 $i$ $(T_i,T_e,T_R)$ $norm(f)$ +-----+-----------------------------------------+------------+ 0 (300.0000000, 300.0000000, 300.0000000) 275.0 1 (288.0000000, 288.0000000, 329.5240078) 287.3 2 (292.4932996, 292.0208497, 314.3715489) 25.47 3 (292.6340330, 292.1467857, 312.7843999) 0.2664 4 (292.6356304, 292.1482151, 312.7673979) 0.00003122 5 (292.6356306, 292.1482153, 312.7673959) 1.137e-12 6 (292.6356306, 292.1482153, 312.7673959) 6.821e-13 7 (292.6356306, 292.1482153, 312.7673959) 6.821e-13 8 (292.6356306, 292.1482153, 312.7673959) 6.821e-13 9 (292.6356306, 292.1482153, 312.7673959) 6.821e-13 10 (292.6356306, 292.1482153, 312.7673959) 6.821e-13 11 (292.6356306, 292.1482153, 312.7673959) 6.821e-13