Python differential analysis of heat loss across a pipe I made a differential equation solver to forecast heat loss through a pipe with heat loss coefficient contours. It does what I want it to do but I had to use global variables to change the outputs of my simulation function - is there a better way to deal with this issue? Another problem I have is that my code is very slow - when I increase dx the answer seems to become less accurate (the graphs seem to approach a single result past 0.0001)

import scipy
# R_w = 0.028# W/mk air heat loss ocefficient
R_w = 0.654# W/mk water heat loss ocefficient
R_f = 0.042 # W/mk pipe heat loss coefficient
r = 0.2 # m pipe radius
T_0 = 50 # C initial temperature
T_a = 20 # C ambiant temperature
dx = 0.0001 # m differential distance
L = 25 # m length of pipe

def heat_balance(T_1,T_2):
   heat_in = (T_1-T_2) * scipy.pi * r**2 * R_w
   heat_out = ( (T_1 + T_2) / 2 - T_a ) * 2 * r * scipy.pi * dx * R_f
   return abs(heat_in-heat_out)

def simulation():
   T = [T_0]
   X = [0]
   print 'r is ',R_f
   import numpy as np
   from scipy.optimize import minimize
   for x in np.arange(dx,L,dx):
     X.append(x)
     fun = lambda s: heat_balance(T[-1],s)
     T.append(float(minimize(fun,T[-1]-dx,method="SLSQP")['x']))
   return([T,X])

from matplotlib import pyplot as plt
typs = ['+','x','|','_','*']
#typs = ['D','s','8','p','^']
cols = ['red','blue','black','purple','brown']

for i, rf in enumerate([7,14,50,100,200]):
    R_f = 1.0/rf
    ans = simulation()
    plt.scatter(ans[1],ans[0],c=cols[i],marker=typs[i],s = 100,label=rf,lw=3)

plt.legend(fontsize=20,title='Pipe Insolation R Value')
plt.xlabel('Distance (m)',fontsize=20)
plt.ylabel('Average Temperature (C)',fontsize=20)
plt.show()
Globals in Python are not necessarily a bad thing, especially for what you're using them for. However, they can be prevented. There are a couple of major flaws I found so I'm glad you came over to get a review. The point you're most worried about will be handled at the end.

Style

Clean Python code adheres to PEP8, the general style guide for Python. Your code violates this guide in a couple of ways. I'll explain the violations I deem most important and fixed it up for you to show the difference. I do recommend you read PEP8 yourself as well.

To verify whether your code adheres to the guide, you can either use the pep8 tool inside Python itself or pep8online. This will validate your code against most of the guide.

Import

Import statements should be on top of your file for maintainability purposes. This is a common practice in many languages, including C, C++, C# and Java. If I now skim your code, I'll see unexpected functions like plt.xlabel() because I didn't notice the following line:

from matplotlib import pyplot as plt


All imports should be on the top.

Globals

Imports should be followed by your global variables. Often those globals are supposed to be constant (note that Python does not actually enforce them being constant). The naming convention for constants is UPPERCASE_WITH_UNDERSCORES. This makes it clear to anyone reading the code the value should not be changed and is probably defined at module level.

This also means that other variables should not use the same casing. All relevant casing styles are in the PEP8.

Variable naming

Now we're handling variables anyway, please use more descriptive variable names. Single-letter variables and functions make for hard-to-read code. Code should be self descriptive. With good names, many of your comments would become redundant. Functions with names are fun can definitely use something more descriptive.

Whitespace

Give your operators and commas some breathing room. foo, bar(woof + 1), fizz is a lot easier on the eyes than foo,bar(woof+1),fizz.

Indentation is done with spaces, not with tabs. Multiples of 4 are preferred, but multiples of 2 are also common.

import scipy
import numpy as np
from scipy.optimize import minimize
from matplotlib import pyplot as plt

R_w = 0.654  # W/mk water heat loss ocefficient
R_f = 0.042  # W/mk pipe heat loss coefficient
PIPE_RADIUS = 0.2  # m pipe radius
PIPE_LENGTH = 25  # m length of pipe
TEMPERATURE_INITIAL = 50  # C initial temperature
TEMPERATURE_AMBIANT = 20  # C ambiant temperature
dx = 0.0001  # m differential distance


def heat_balance(T_1, T_2):
    heat_in = (T_1-T_2) * scipy.pi * PIPE_RADIUS**2 * R_w
    heat_out = ( (T_1 + T_2) / 2 - TEMPERATURE_AMBIANT ) * 2 * PIPE_RADIUS * scipy.pi * dx * R_f
    return abs(heat_in-heat_out)


def simulation():
    T = [TEMPERATURE_INITIAL]
    X = [0]
    print 'r is ', R_f
    for x in np.arange(dx, PIPE_LENGTH, dx):
        X.append(x)
        fun = lambda s: heat_balance(T[-1], s)
        T.append(float(minimize(fun, T[-1] - dx, method="SLSQP")['x']))
    return([T, X])

typs = ['+', 'x', '|', '_', '*']
cols = ['red', 'blue', 'black', 'purple', 'brown']

for i, rf in enumerate([7, 14, 50, 100, 200]):
    R_f = 1.0/rf
    ans = simulation()
    plt.scatter(ans[1], ans[0], c=cols[i], marker=typs[i], s = 100, label=rf, lw=3)

plt.legend(fontsize=20, title='Pipe Insolation R Value')
plt.xlabel('Distance (m)', fontsize=20)
plt.ylabel('Average Temperature (C)', fontsize=20)
plt.show()


Objective Oriented Programming

This is one of the best solutions against global variables: objective oriented programming.

Did you know you can make classes in Python?

class PipeCalculator:
    self.R_w = 0.654  # W/mk water heat loss ocefficient
    self.R_f = 0.042  # W/mk pipe heat loss coefficient
    self.PIPE_RADIUS = 0.2  # m pipe radius
    self.PIPE_LENGTH = 25  # m length of pipe
    self.TEMPERATURE_INITIAL = 50  # C initial temperature
    self.TEMPERATURE_AMBIANT = 20  # C ambiant temperature
    self.dx = 0.0001  # m differential distance


    def heat_balance(self, T_1, T_2):
        heat_in = (T_1-T_2) * scipy.pi * self.PIPE_RADIUS**2 * R_w
        heat_out = ( (T_1 + T_2) / 2 - self.TEMPERATURE_AMBIANT ) * 2 * self.PIPE_RADIUS * scipy.pi * dx * R_f
        return abs(heat_in-heat_out)


    def simulation(self):
        T = [self.TEMPERATURE_INITIAL]
        X = [0]
        print 'r is ', self.R_f
        for x in np.arange(self.dx, self.PIPE_LENGTH, self.dx):
            X.append(x)
            fun = lambda s: self.heat_balance(T[-1], s)
            T.append(float(minimize(fun, T[-1] - self.dx, method="SLSQP")['x']))
        return([T, X])

    def plot_calculations(self):
        typs = ['+', 'x', '|', '_', '*']
        cols = ['red', 'blue', 'black', 'purple', 'brown']

        for i, rf in enumerate([7, 14, 50, 100, 200]):
            self.R_f = 1.0/rf
            ans = simulation()
            plt.scatter(ans[1], ans[0], c=cols[i], marker=typs[i], s = 100, label=rf, lw=3)

        plt.legend(fontsize=20, title='Pipe Insolation R Value')
        plt.xlabel('Distance (m)', fontsize=20)
        plt.ylabel('Average Temperature (C)', fontsize=20)
        plt.show()

calculator = PipeCalculator()
calculator.plot_calculations()


And poof, all the globals are gone.

Of-course, a straight port from what you have to object oriented still looks ugly. Pass some variables around, use the power of arguments and when in doubt don't hesitate to split up to more functions. But at least you now have a start and it will be easier to re-use this code later on in different projects.

As two final notes:


Your lambda could be re-written to a normal function for improved readability. This may help you figure out how to optimize the calculations.
Pipes usually have an inside and outside diameter. You may want to explicitly state which one you're using.