I am trying to carry out an iterative calculation in python using a finite difference method. I found the finite difference method from this:

http://depa.fquim.unam.mx/amyd/archivero/DiferenciasFinitas3_25332.pdf

The values that the code calculates are correct. The problem is that it is only showing the final values. What I want is to extract the values for any point in the x-direction so that I can plot them, as well as extract the values at any point in time e.g the values of the points halfway through the calculations. Is this the right way for doing an iterative calculation? The code is shown below:

import numpy as np
import scipy as sp
import time
import matplotlib as p
L=0.005
Nx=3
T=5
N1=5
k=0.5
rho=1200
c=1000
a=(k/(rho*c))
x = np.linspace(0, L, Nx+1)   # mesh points in space
dx = x[1] - x[0]
t = np.linspace(0, T, N1)    # time
dt = t[1] - t[0]
toutside=5
Coefficient = a*dt/dx**2
bi=0.5
ui = sp.zeros(Nx+1)
u = sp.zeros(Nx+1)
for i in range(Nx+1):
  ui[i] = 50 # initial values

for n in range(0, N1):

   for i in range(0,1):
      u[i] = 2*Coefficient*(ui[i+1]+bi*toutside)+(1-2*Coefficient-2*bi*Coefficient)*ui[i]
   for i in range(1,Nx):
      u[i] = Coefficient*(ui[i+1]+ui[i-1])+(1-2*Coefficient)*ui[i]
   for i in range(Nx,Nx+1):
      u[i] = 2*Coefficient*(ui[i-1])+(1-2*Coefficient)*ui[i]
      ui[:]= u #updates matrix for next loop
 print ui

I've modified my code to based on danodonovan answer to:

for n in range(0, N1):

for i in range(0,1):
    u[i] = 2*Coefficient*(ui[i+1]+bi*toutside)+(1-2*Coefficient-2*bi*Coefficient)*ui[i]
for i in range(1,Nx):
    u[i] = Coefficient*(ui[i+1]+ui[i-1])+(1-2*Coefficient)*ui[i]
for i in range(Nx,Nx+1):
    u[i] = 2*Coefficient*(ui[i-1])+(1-2*Coefficient)*ui[i]
    ui=u
    a=list(ui)

print a

When I try to take the entire list out of the loop, only the final values are produced. How do I extract the whole list? Is this the right way to do an iterative calculation using the values of the previous row to calculate the values of the new row?

有帮助吗?

解决方案

For a starter indent the print ui to see the values of ui for all N1. Then append your results to a list:

res = []
for n in range(0, N1):

   for i in range(0,1):
      u[i] = 2*Coefficient*(ui[i+1]+bi*toutside)+(1-2*Coefficient-2*bi*Coefficient)*ui[i]
   for i in range(1,Nx):
      u[i] = Coefficient*(ui[i+1]+ui[i-1])+(1-2*Coefficient)*ui[i]
   for i in range(Nx,Nx+1):
      u[i] = 2*Coefficient*(ui[i-1])+(1-2*Coefficient)*ui[i]
      ui[:]= u #updates matrix for next loop
   print ui
   res.append(ui.copy())

print res

Produces this result:

[array([ 41.5625,  50.    ,  50.    ,  50.    ]),
 array([ 37.87109375,  48.41796875,  50.        ,  50.        ]),
 array([ 35.6628418 ,  46.73706055,  49.70336914,  50.        ]),
 array([ 34.06639099,  45.21682739,  49.20280457,  49.88876343]),
 array([ 32.79785633,  43.87349129,  48.58405113,  49.63152885])]

其他提示

(If I understand your question correctly) Using matplotlib you can do

 import matplotlib as p

 # and after your loop has completed

 p.pyplot.plot(range(0, N1), ui, 'o-')
 p.pyplot.show()

to get a simple plot of your data u against range(0, N1).

I don't know what you expect ui to be ui[:]= u is an odd thing to be doing - it sets a copy of ui to u but you don't retain the copy of ui that ui[:] is producing.

Hint: think about ui[:] as being the same as list(ui)

许可以下: CC-BY-SA归因
不隶属于 StackOverflow
scroll top