Help with Homework

[rjmiller]

Ok, So i'm attempting to simulate diffusion through a 2-D plate of multiple concentrations of compounds and my program so far looks something like follows

from scipy import *
from pylab import *

#==============================================================================
#simulation conditions
#==============================================================================
nx = 64                     # The number of x nodes
ny = 64                     # The number of y nodes
h = 0.01                    #spacing between nodes
Du = 2*10**-5               # Diffusion coefficient of U
Dv = 10**-5                 # Diffusion coefficient of V
delta_t = 0.9*h**2/(4*Du)   # time step


def lap_2d(f, h):
    '''Return the laplacian of a 2-D array f with node spacing h.'''
    return (roll(f, -1, axis=0) + roll(f, 1, axis=0) - 4*f + roll(f, -1, axis=1) + roll(f, 1, axis=1))/(h**2)

def dstatedt(F, k, U, V):
    return ((Du*f(U, h)) - U*V**2 + F*(1-U) & ((Dv*f(V, h)) + U*V**2 - (k + F)*V)

#==============================================================================
# Initial conditions
#==============================================================================
U = zeros((nx, ny), dtype = float) # a 2-D array 
V = zeros((nx, ny), dtype = float) # a 2-d array 

U[(x > 12*nx/16) & (x < 6*nx/16)] = 1.0 +0.01* stats.uniform.rvs() # original U concentration 
V[ (y > 12*ny/16) & (y < 6*ny/16)] = 0.0 + 0.01*stats.uniform.rvs() # original V concentration 

U[(x > 6*nx/16 ) & (x < 12*nx/16)] = 0.5 + 0.01*stats.uniform.rvs() # initial U concentration added
V[(y > 6*ny/16 ) & (y < 12*ny/16)] = 0.25 + 0.01*stats.uniform.rvs() #initial V concentration added

If you notice I am running into several issues. With my function dstatedt I want it to return two different numbers in some sort of an array, because where as both are based on concentrations of U and V, they are changing at different rates. However, I want them written as one function.

also, the last lines of code talk about the initial and original concentrations. As you can see I want to added random noise type data to the arrays, but only for certain parts of the arrays. My teachers gave me hints to use the stats.uniform.rvs() function, but i don't know how to manipulate the loc, scale, and size arguments so that i can get +- 0.01 random noise.

If anyone has any sort of clues as to how to fix my issues, that would be great.

[crono5788]

I wish I knew the subject matter, I might be able to help you better. I'll try anyway, though :)

For your first issue, you can return a tuple or a list:
def someFunction:
return ("Value1", "Value2")
Or, alternatively, for a list (which is mutable, unlike a tuple)
return ["Value1", "Value2"]

And in your code...
x = someFuction()
print x[0] # Value1
print x[1] # Value2

For your second issue, do you need only the values 0.01 or -0.01, or do you want a value x, so that -0.01 <= x <= 0.01?

[ithelp]

I guess there is few floating point arithmetic erros as well.

[rjmiller]

Alright, I figured it out, I'll post my code below if any of you care to take a look:

from scipy import *
from pylab import *

#==============================================================================
# simulation conditions
#==============================================================================
nx = 64                     # The number of x nodes
ny = 64                     # The number of y nodes
h = 0.01                    # spacing between nodes
Du = 2e-5                   # Diffusion coefficient of U
Dv = 1e-5                   # Diffusion coefficient of V
delta_t = 0.9*h**2/(4*Du)   # time step
time_step = 10000           # time iterations

#==============================================================================
# Initial conditions
#==============================================================================
U = ones((nx, ny), dtype = float)   # a 2-D array 
V = zeros((nx, ny), dtype = float)  # a 2-d array 


U[24:40, 24:40] = 0.5 
V[24:40, 24:40] = 0.25

U += stats.uniform.rvs(size = (nx, ny), loc = -.005, scale = 0.01) # original U concentration 
V += stats.uniform.rvs(size = (nx, ny), loc = -.005, scale = 0.01) # original V concentration
 

#==============================================================================
# Unique k and F values
#==============================================================================

##k = 0.060
##F = 0.040

##k = 0.056
##F = 0.024

k = 0.059
F = 0.040

##k = 0.056
##F = 0.020

##k = 0.050
##F = 0.016


def lap_2d(f, h):
    '''Return the laplacian of a 2-D array f with node spacing h.'''
    return (roll(f, -1, axis=0) + roll(f, 1, axis=0) - 4*f + roll(f, -1, axis=1) + roll(f, 1, axis=1))/(h**2)

def dstatedt(U, V, F, k):
    '''Returns a tuple of arrays with new data interpreted given values of F and k'''
    uvsq = U*V**2
    du = Du*lap_2d(U, h) - uvsq + F*(1-U)
    dv = Dv*lap_2d(V, h) + uvsq - (k + F)*V
    return (du, dv)

for step in range(time_step):
    du, dv = dstatedt(U, V, F, k)
    U += delta_t*du
    V += delta_t*dv
    if step % 10 == 0: # saves image every 10 time steps 
        subplot(121)
        imshow(U, vmin=0, vmax = 1)
        subplot(122)
        imshow(V, vmin = 0, vmax = 1)
        filename = 'image%06d.png' % step
        clf()
from os import *

system(r"c:\mplayer\mencoder.exe -ovc lavc -lavcopts vcodec=xvid:vbitrate=500 -mf type=png:fps=10 -nosound mf://*.png -o test.avi")