How do I implement a 1D linear convection using finite differences in Python?

In summary, the programmer is trying to implement a simple one dimensional linear convection problem in Python, but is struggling with understanding how to get started. They have Python 2.6, Numpy and matplotlib installed, and have attached the pseudocode as an image file. They are not sure how to start the loops with a value of 'i' and maybe need to add and increment to 'i' once it has done the loop once. They are also unsure of what the nested loop is supposed to be doing.
  • #1
Jonny6001
20
0
Hi, I have a pseudocode I would like to try and implement and understand the programming of it in python. The physics and maths of the case is no problem, but the implementation of the code in Python is something I'm not familiar with.

The problem is a simple 1D linear convection using finite differences.

I have Python 2.6, Numpy and matplotlib.

I have attached the pseudocode as an image file. I can more or less input this in a form that Python understands but I'm not quite sure how to start off the loops with a value of 'i' and maybe I need to add and increment to 'i' once it has done the loop once?
I mean for a start, I should be able to run the first for loop which says for values of i between 0.5 and 1, the value of 'ui' will be 2, everywhere else ui=1.
Then plot the results in matplotlib.

I have tried this and I get a blank ploot with x-axis around the 20 and y-axis around 1, so the axes are not far off but the increment of the loop is something I'm not getting. By the way, this is my own interest, not a homework or assignment, so there isn't anyone I can ask the most basic quiestions to.

Thanks a lot for your time, I'm sure if I can be shown how to make a start then it will slowly click into place.
 

Attachments

  • Untitled.jpg
    Untitled.jpg
    10.2 KB · Views: 559
Technology news on Phys.org
  • #2
Can you please enter your code in your post, not as an image? I can either look at your code or type a response, but not both.
 
  • #3
nx=20, nt=50
dt=0.01, c=1
dx=2/(nx-1)

for i in (1,nx):
if 0.5<=x(i)<=1:
u(i)=2
else:
u(i)=1
end
end

for it in (1,nt):
un=u
for i in (2,nx-1):
u(i)=un(i)-c*dt/dx(un(i)-un(i-1))
end
end


Sorry about that, it didn't cross my mind.

Thank you
 
  • #4
Jonny6001 said:
Formatted your code using [ code] [ /code] block (without the leading space).
Code:
nx=20, nt=50
dt=0.01, c=1
dx=2/(nx-1)

for i in (1,nx):
	if 0.5<=x(i)<=1: // see comment 1 below
		u(i)=2
	else:
		u(i)=1
	end
end

for it in (1,nt):
	un=u  // see comment 2 below
	for i in (2,nx-1):
		u(i)=un(i)-c*dt/dx(un(i)-un(i-1)) // see comment 3
	end
end
[\code]

Sorry about that, it didn't cross my mind.

Thank you[/QUOTE]

Caveat: I don't know python, but I'm very familiar with C, C++, C#, and a slew of other programming languages.

1. I'm 99% sure you can't do this.  0.5<=x(i)<=1 is strictly mathematics notation shorthand. Your intent is clear, but you can't do it this way. To check that x(i) is between 0.5 and 1.0, do this: 
if ( x(i) >= 0.5 && x(i) <= 1.0) 
{
     ...
}
else
{
    ...
}
What I've written is C/C++ syntax, so you'll need to modify it a bit for Python. The else clause is executed when the expression in the if clause is false; i.e., when x(i) < 0.5 OR when x(i) > 1.0. If you need to take different actions for these then you need an elseif (or maybe it's elif in python - don't know) and an else clause.

2. What is un = u doing? Both of these appear to be arrays. If arrays in python work anything like they do in C, what this is doing is to put the address of the u array into un.

3. In this line, u(i)=un(i)-c*dt/dx(un(i)-un(i-1)), you're missing a multiplication operator - *. It should be u(i)=un(i)-c*dt/dx * (un(i)-un(i-1)) 

[quote=Jonny6001] I can more or less input this in a form that Python understands but I'm not quite sure how to start off the loops with a value of 'i' and maybe I need to add and increment to 'i' once it has done the loop once?
I mean for a start, I should be able to run the first for loop which says for values of i between 0.5 and 1, the value of 'ui' will be 2, everywhere else ui=1.
[/quote]
I don't understand what you're asking about the loop variables. In your first loop i takes on values of 1, 2, 3, ..., 19, 20. In the second loop in takes on values of 1, 2, 3, ..., 50. 

Also, in your first loop, if [B]x(i)[/B] (not i) is between .5 and 1, [B]u(i)[/B] (not ui) gets set to 2; otherwise [B]u(i)[/B] gets set to 1. Try to be more precise in what you say.

Finally, in your nested loop, what are you trying to do? The outer loop runs 50 times, and each iteration of the outer loop causes the inner loop to run 18 times. This means that each statement in the inner loop runs 50 * 18 = 900 times.
 
  • #5
Hi, thanks for taking the time to reply. I am pretty sure that in python, you can do 0.5<=x(i)<=1, I have seen a few documents and it doesn't complain about the syntax.

The issue I think is the loop doesn't seem to be running 'i' between the values in the brackets.

The code below, I see it as being; increment 'i' between 1 and 'nx', if i<10, put y=10, if y==10 put y=20, if y>10 put y=50

So if I was to plot (i,y) I was expecting to see a graph with 'i' along the x-axis ranging from 1-20, and the corresponding values of 'y' so 'y' would be 10, then when i==10, y goes up to 20 then its 50 for the remaining values of 'i'.

For some reason the plot axis goes to the region of the final value as it it wasn't running and recording the whole loop.

I won't paste the code in this post as its so short, you can see the code and the graph output in the attached image.

I recall ages ago when I did more things like this, there was a commond at the end of the loop such as i++ or something to add 1 to the value of i, perhaps Python requires something similar?

I'm sorry if my questions are a little boring but I really want to get heavily into Scientific Programming and Mathematical Modelling.
By the way, the pseudocode I wrote earlier is a way of approximating a partial differential equation describing linear convection, using finite difference method.

Thanks so much for taking the time to bother with this.
 

Attachments

  • Untitled1.jpg
    Untitled1.jpg
    25.5 KB · Views: 473
  • #6
Jonny6001 said:
Hi, thanks for taking the time to reply. I am pretty sure that in python, you can do 0.5<=x(i)<=1, I have seen a few documents and it doesn't complain about the syntax.
Nor does C, which is even worse. IOW, it's not a syntax error in C, so you don't hear from the compiler about it, but it is a semantic error.

Here's an example that shows how it works in C and might work in python, using a slightly simpler example.

result = 0.5 <= x <= 1

Suppose x has been set to .75. The expression 0.5 <= x is true, so the value of the expression (0.5 <= x) is 1. 1 <= 1 is true, so the overall expression has a value of 1, which is what is stored in result. That is the expected behavior.

Now suppose x has been set to 1.5. The expression 0.5 <= x is true again, so the value of the expression (0.5 <= x) is 1. 1 <= 1 is true, so the overall expression has a value of 1, which is what is stored in result. This is not what was expected.

One more. Suppose x has been set to 0.0. The expression 0.5 <= x is false this time, so the value of the expression (0.5 <= x) is 0. 0 <= 1 is true, so the overall expression has a value of 1, which is what is stored in result. This is not what was expected.

In other words, no matter what value x is, the expression 0.5 <= x <= 1 always comes out true - at least in C, and I would be very surprised if python behaves any differently.

I'll take a look at your other comments in a little while.

Jonny6001 said:
The issue I think is the loop doesn't seem to be running 'i' between the values in the brackets.

The code below, I see it as being; increment 'i' between 1 and 'nx', if i<10, put y=10, if y==10 put y=20, if y>10 put y=50

So if I was to plot (i,y) I was expecting to see a graph with 'i' along the x-axis ranging from 1-20, and the corresponding values of 'y' so 'y' would be 10, then when i==10, y goes up to 20 then its 50 for the remaining values of 'i'.

For some reason the plot axis goes to the region of the final value as it it wasn't running and recording the whole loop.

I won't paste the code in this post as its so short, you can see the code and the graph output in the attached image.

I recall ages ago when I did more things like this, there was a commond at the end of the loop such as i++ or something to add 1 to the value of i, perhaps Python requires something similar?

I'm sorry if my questions are a little boring but I really want to get heavily into Scientific Programming and Mathematical Modelling.
By the way, the pseudocode I wrote earlier is a way of approximating a partial differential equation describing linear convection, using finite difference method.

Thanks so much for taking the time to bother with this.
 
  • #7
Jonny6001 said:
Hi, thanks for taking the time to reply. I am pretty sure that in python, you can do 0.5<=x(i)<=1, I have seen a few documents and it doesn't complain about the syntax.

The issue I think is the loop doesn't seem to be running 'i' between the values in the brackets.

The code below, I see it as being; increment 'i' between 1 and 'nx', if i<10, put y=10, if y==10 put y=20, if y>10 put y=50

So if I was to plot (i,y) I was expecting to see a graph with 'i' along the x-axis ranging from 1-20, and the corresponding values of 'y' so 'y' would be 10, then when i==10, y goes up to 20 then its 50 for the remaining values of 'i'.
For matlab's plot function, the two arguments have to be arrays, and in your code they don't seem to be arrays.

Here's a python tutorial that I found. Take a look at section 4.2, which talks about for loops. You don't need to increment your loop counter.

Also, I looked at the section on logical conditions and didn't see anything to contradict what I said in my previous reply about what 0.5 <= x(i) <= 2 would evaluate as. I still believe this always evaluates to 1 (or true).
http://docs.python.org/tutorial/
Jonny6001 said:
For some reason the plot axis goes to the region of the final value as it it wasn't running and recording the whole loop.

I won't paste the code in this post as its so short, you can see the code and the graph output in the attached image.

I recall ages ago when I did more things like this, there was a commond at the end of the loop such as i++ or something to add 1 to the value of i, perhaps Python requires something similar?

I'm sorry if my questions are a little boring but I really want to get heavily into Scientific Programming and Mathematical Modelling.
By the way, the pseudocode I wrote earlier is a way of approximating a partial differential equation describing linear convection, using finite difference method.

Thanks so much for taking the time to bother with this.
 
  • #8
Mark44 said:
In other words, no matter what value x is, the expression 0.5 <= x <= 1 always comes out true - at least in C, and I would be very surprised if python behaves any differently.

Python is a very different language to C! Although Python (the standard version, at least) is written in C, the style, intent, and scope of the language is extremely unlike what you might be used to in C.

For what it's worth, the change in mindset that I needed to undergo when I first learned Python was quite painful; I had a solid C/C++ background and many Pythonic ways of doing things just seemed wrong. With the benefit of hindsight though, it's clear that Python gives rise to a much more natural style of expression than C/C++.

Mark44 said:
Also, I looked at the section on logical conditions and didn't see anything to contradict what I said in my previous reply about what 0.5 <= x(i) <= 2 would evaluate as. I still believe this always evaluates to 1 (or true).

You can have several arithmetic tests as conditions to an if statement in Python. For instance, consider the following situation where I define a Python list which contains the integers from one to ten, inclusive.

Code:
# Define a list containing the integers from 1 to 10.
li = range(1, 11)

# Print out whether the numbers in the list are between 6 and 9, inclusive.
for x in li:
    if 6 <= x <= 9:
        print x
    else:
        print 'FAIL!'

The output from this is as one would expect:

Code:
FAIL!
FAIL!
FAIL!
FAIL!
FAIL!
6
7
8
9
FAIL!
 
Last edited:
  • #10
Ok, I have tried to take a step back to better understand, it's not enough to be good at the mathematics I'm trying to express in a program, I understand that now.

I have attached something I wrote to see how the loops work and how the plot function in matplotlib works.
It appears that the for loop runs ok and increments up to value of 'nx', and as you can see I have an array of 'i' and corresponding 'y' values. I thought that if I plotted (i,y) this should plot the whole range, but it seems it sort of plots the last value of (i,y) in the array?

Does anyone know how I should plot the values, so between 5<=i<=10, the value of 'y' should be 2, everywhere else it's 1.

Thanks a lot.
 

Attachments

  • Untitled3.jpg
    Untitled3.jpg
    23.6 KB · Views: 447
  • #11
Jonny6001 said:
Ok, I have tried to take a step back to better understand, it's not enough to be good at the mathematics I'm trying to express in a program, I understand that now.

I have attached something I wrote to see how the loops work and how the plot function in matplotlib works.
It appears that the for loop runs ok and increments up to value of 'nx', and as you can see I have an array of 'i' and corresponding 'y' values. I thought that if I plotted (i,y) this should plot the whole range, but it seems it sort of plots the last value of (i,y) in the array?

Does anyone know how I should plot the values, so between 5<=i<=10, the value of 'y' should be 2, everywhere else it's 1.

Thanks a lot.


Instead of attaching your code as a jpeg, it's much better to use the [ code ] [ /code ] tags (although without the spaces inside the brackets!). For instance, the code you've given in your last attachment is as follows.

Code:
import matplotlib.pyplot as plt

nx = range(1, 21)

for i in nx:
	if 5 <= i <= 10:
		
		y = 2
		print i, y
	
	else: 
		
		y = 1
		print i, y


plt.plot(i, y)
plt.show()

The intent behind what you're doing is almost correct, but the implementation is slightly wrong. In each step of the for loop you run an if..else test. The condition for the test is that [itex]5 \le i \le 10[/itex]; if this is TRUE, you define [itex]y=2[/itex] and print the values of [itex]i[/itex] and [itex]y[/itex] to standard output. You then go on to the next stage in the loop and run the test again. The important thing to note is that when you go to the next step in the loop you are (i) not storing the values of [itex]i[/itex] and [itex]y[/itex] you had in the previous step, and (ii) you are rewriting the old values with new ones.

What you really want to do is to store the appropriate values somewhere so that they can be plotted later on. Here's a simple example of how this can be achieved.

Code:
import matplotlib.pyplot as plt

nx = range(1, 21)

# Define lists to hold the values of i and y
iList = []
yList = []

# Now run through the loop.
for i in nx:
	
	# Add the current value of i to the iList
	iList.append(i)
	
	if 5 <= i <= 10:
		
		# Append 2 to the yList
		yList.append(2)
	
	else: 
		
		# Append 1 to the yList
		yList.append(1)


# Now print out the contents of the lists as ordered pairs.
for elem in range(len(iList)):
	print '(%g, %g)' % (iList[elem], yList[elem])


# Now plot the data contained in the lists. Note that you can pass other options
# to the plot command to give you, for example, a plot of just the data points
# with no interpolating line between them. The relevant command for this would
# be something of the form:
# plt.plot(iList, yList, 'bo') 
plt.plot(iList, yList)
plt.show()

This results in the pairs of numbers being printed to standard output and will give you the following plot.

[PLAIN]http://img686.imageshack.us/img686/8480/plotx.png [Broken]

It's possible to do this more elegantly, but it's good to concentrate on getting code working before you attempt to make it beautiful.
 
Last edited by a moderator:
  • #12
Brilliant, thanks a lot for the reply, I clearly see how this is achieved now.

I will try and get the methodology working in a more practical way, to try and plot the values of the initial case.
Have you any experience in finite difference methods in Python, such as Burgers Equation, Convection equation or what I eventually want to build up to, Navier Stokes.

Again, thanks for your time.
 
  • #13
Jonny6001 said:
Have you any experience in finite difference methods in Python, such as Burgers Equation, Convection equation or what I eventually want to build up to, Navier Stokes.

Unfortunately, I know very little about finite difference schemes. However, a quick google reveals plenty of people who are doing just this sort of thing. You can find course notes on numerical methods in Python, including a code for a finite differences approach to the one-dimensional convection equation http://www.maths.uq.edu.au/courses/MATH3203/ [Broken], for example. I also imagine plenty of others here will have lots of experience with this and will be able to answer any questions you might have.
 
Last edited by a moderator:

1. What is 1D Linear Convection in Python?

1D Linear Convection in Python refers to the simulation of the transport of a physical quantity, such as heat or fluid flow, in a single dimension using the Python programming language. It involves solving a set of equations that describe the behavior of the quantity being transported in a linear manner.

2. What is the difference between 1D Linear Convection and 2D/3D Convection?

The main difference between 1D Linear Convection and 2D/3D Convection is the number of dimensions in which the physical quantity is being transported. 1D Linear Convection only considers transport in a single dimension, while 2D/3D Convection takes into account transport in two or three dimensions.

3. How is 1D Linear Convection simulated in Python?

1D Linear Convection can be simulated in Python using a numerical method known as Finite Difference Method (FDM). This involves discretizing the domain into a grid and approximating the differential equations that describe the behavior of the physical quantity being transported using finite differences.

4. What are the applications of 1D Linear Convection in Python?

1D Linear Convection in Python has various applications in the fields of fluid mechanics, heat transfer, and atmospheric sciences. It can be used to model the transport of pollutants in the atmosphere, the flow of heat in materials, and the movement of fluids in pipes or channels.

5. Are there any limitations to using 1D Linear Convection in Python?

One limitation of using 1D Linear Convection in Python is that it only considers transport in a single dimension, which may not accurately represent real-world scenarios. Additionally, the numerical method used for simulation may introduce errors and may not always converge to an accurate solution.

Similar threads

  • Programming and Computer Science
Replies
1
Views
874
  • Programming and Computer Science
Replies
4
Views
4K
  • Programming and Computer Science
Replies
2
Views
1K
Replies
6
Views
2K
  • Programming and Computer Science
Replies
2
Views
2K
  • Programming and Computer Science
Replies
3
Views
2K
  • Programming and Computer Science
Replies
2
Views
16K
  • Programming and Computer Science
Replies
5
Views
4K
  • Programming and Computer Science
Replies
4
Views
6K
Back
Top