How can I limit the array in MathCAD to prevent crashing?

[haitrungdo82]

Homework Statement

I don't know how to limit the array. I know that the program crashed because it exceeds the array end points at zero and maximum position and time. But how can I tell MathCAD not to do so? I'm stuck without any idea! Can anyone who can show me the code lines I should include. Thank you very much!

Homework Equations

The Attempt at a Solution

The following is what I have achieved so far. I would really appreciate any input from you! Thank you very much!

Δt := 1 min
Δx := 5 m

v := 16.09 km/h
D := 0.17*10^-4 m^2/min

t := 0 ... 1999
x := 1 ... 3999
f_t,x := 1

f_t+1,x := f_t,x - (v*Δt/2)*(f_t,x+1 - f_t,x-1) + [D*Δt/(Δx)^2]*(f_t,x-1 - 2*f_t,x

+ f_t,x+1)

Then it crashed! Can anyone help me, please?

 

[gneill]

t := 0 ... 1999
x := 1 ... 3999
f_t,x := 1

f_t+1,x := f_t,x - (v*Δt/2)*(f_t,x+1 - f_t,x-1) + [D*Δt/(Δx)^2]*(f_t,x-1 - 2*f_t,x

+ f_t,x+1)

Then it crashed! Can anyone help me, please?

It looks like you're trying to evaluate 8 million array items. It your machine has enough memory, it probably didn't crash but rather the mathcad application has gone catatonic while it grinds away at the numbers (for each t value it's evaluating 4000 x's).

Reduce your t and x ranges to something trivial for testing purposes and check to see if what its doing is really what you expect it to be doing.

 

[haitrungdo82]

Thank you for your input, qneil!

I reduced the range of x and t to something like t := 199 and x := 299, but the problem doesn't go away.

When I ran the program or managed to correct the mistakes, either "x+1" or "t+1" indices were highlighted in red. That's why I think I exceeded the limit. But I have no idea how to handle this kind of situation! Is there any suggestion?

 

[haitrungdo82]

It looks like you're trying to evaluate 8 million array items. It your machine has enough memory, it probably didn't crash but rather the mathcad application has gone catatonic while it grinds away at the numbers (for each t value it's evaluating 4000 x's).

Reduce your t and x ranges to something trivial for testing purposes and check to see if what its doing is really what you expect it to be doing.

Thank you for your input, qneil!

I reduced the range of x and t to something like t := 199 and x := 299, but the problem doesn't go away.

When I ran the program or managed to correct the mistakes, either "x+1" or "t+1" indices were highlighted in red. That's why I think I exceeded the limit. But I have no idea how to handle this kind of situation! Is there any suggestion?

 

[gneill]

Δt := 1 min
Δx := 5 m

v := 16.09 km/h
D := 0.17*10^-4 m^2/min

t := 0 ... 1999
x := 1 ... 3999
f_t,x := 1

f_t+1,x := f_t,x - (v*Δt/2)*(f_t,x+1 - f_t,x-1) + [D*Δt/(Δx)^2]*(f_t,x-1 - 2*f_t,x

+ f_t,x+1)

First, your term (v*Δt/2) is going to have units of meters, while the rest of the terms would appear to be unit-less. This is going to result in a units mismatch error. Perhaps this rendering of the program is not correct?

You are also pre-filling the entire f array with 1's, then you proceed to overwrite them all but the first one, f_0,1. Any particular reason for performing all the extra work?

You are using an index of x-1. Since x is defined to go from 1..3999, it will reference f_0.0 which does not exist. You will have to set a starting value there.

When you say that MathCAD crashed, do you mean it locked up (became unresponsive or disappeared altogether), or you you mean that one or more expressions are reporting errors?

 

[haitrungdo82]

Δt := 1 min
Δx := 5 m

v := 16.09 km/h
D := 0.17*10^-4 m^2/min

t := 0 ... 1999
x := 1 ... 3999
f_t,x := 1

f_t+1,x := f_t,x - (v*Δt/2)*(f_t,x+1 - f_t,x-1) + [D*Δt/(Δx)^2]*(f_t,x-1 - 2*f_t,x

+ f_t,x+1)

First, your term (v*Δt/2) is going to have units of meters, while the rest of the terms would appear to be unit-less. This is going to result in a units mismatch error. Perhaps this rendering of the program is not correct?

You are also pre-filling the entire f array with 1's, then you proceed to overwrite them all but the first one, f_0,1. Any particular reason for performing all the extra work?

You are using an index of x-1. Since x is defined to go from 1..3999, it will reference f_0.0 which does not exist. You will have to set a starting value there.

When you say that MathCAD crashed, do you mean it locked up (became unresponsive or disappeared altogether), or you you mean that one or more expressions are reporting errors?

Hi qneill,

You made really good comments!

This programming is part of a bigger question. The question wants me to simulate the dipersion of airborne pollutants. More specifically, I have to model numerically the time dependent spatial dispersion of the one-dimensional probability density f(t,x). v is the velocity in x direction of the wind that carries the polluting gas, D is the difussion coefficient. The program requires that, for any given time t, f_x is calculated for all distances x. Then t is increased by D_t, and again f_x is calculated for all distances x, and so on.

f(t=0,x=x_0) = 1 is the initial condition. You're right! I'm also very confused about this initial condition. I think it should be f_t,x-1 = 1, right? But if you look at the attachment, which was copied directly from my Professors's solution which he gave me as a hint, he put f_t,x. May be I misunderstood or misidentified the initial point. That's why I'm very confused!

I noticed the unit inconsistence. But first, I want to follow my Professor's hint (the attachment), then I will go back and take care of the units later. The attached equation is the numerical version of the diffusion equation ( or Fokker-Planck equation).

I think we arrive at the same point of confusion now! Yes, I was also wondering how to define the end points. I think in order to limit the array, all the end points must be defined. Neither the question nor my Professor's hint mentioned about the conditions of the end points, except the initial condition f(t=0,x=x_0) = 1. I think they are implied in the question, but I just can not get it! If I can successfully define the end points, I think I can solve the rest of the question.

When I said "crash", I meant the expression was highlighted in red, i.e. the x+1 or the t+1 indices was highlighted in red. I'm so confused that I don't know which ones I miss and which ones are incorrect.

If you would like more information, I would be happy to give you more. I was just afraid that my post will be deleted if I post the whole question.

 

[haitrungdo82]

Δt := 1 min
Δx := 5 m

v := 16.09 km/h
D := 0.17*10^-4 m^2/min

t := 0 ... 1999
x := 1 ... 3999
f_t,x := 1

f_t+1,x := f_t,x - (v*Δt/2)*(f_t,x+1 - f_t,x-1) + [D*Δt/(Δx)^2]*(f_t,x-1 - 2*f_t,x

+ f_t,x+1)

First, your term (v*Δt/2) is going to have units of meters, while the rest of the terms would appear to be unit-less. This is going to result in a units mismatch error. Perhaps this rendering of the program is not correct?

You are also pre-filling the entire f array with 1's, then you proceed to overwrite them all but the first one, f_0,1. Any particular reason for performing all the extra work?

You are using an index of x-1. Since x is defined to go from 1..3999, it will reference f_0.0 which does not exist. You will have to set a starting value there.

When you say that MathCAD crashed, do you mean it locked up (became unresponsive or disappeared altogether), or you you mean that one or more expressions are reporting errors?

Hi qneill,

Btw, if you would like to email me, feel free to do so. My email is << e-mail address deleted by Moderators >>Thank you!

 

[gneill]

Suppose you define your index variables as:

nt := 2000
nx := 4000

t := 0 .. nt - 1
x = 1 .. nx - 1

Your difference equation references f_t+1,x on the LHS, while on the RHS it needs values for f_t,x , f_t,x-1 , f_t,x+1 . This means that it is looking for values calculated at the previous timestep, t = 0, when it's assigning values to f_t+1,__. That's fine.

Now check the x indexes. On the LHS they range from 1 though 3999. On the right they want to go from 0 (for f_t,x-1) through 4000 (for f_t,x+1) at the ends of the range. That means you're going to need to pre-assign values for x=0 and x=4000 for all of t.

Presumably, since you mention that this is a diffusion problem, these pre-assigned values will represent the boundary values for the diffusion field.

I would suggest, then, setting up the entire first row of your f_0,x with the initial concentration pattern, and then setting:

Cb = {some value representing the concentration at the boundaries (could be zero)}
f_t,0 = Cb
f_t,nx = Cb

and then proceed with the diffusion calculations.

You're going to have to work out the issues with the weighting term v*Dt/2. I don't think it should have units or be greater than 1.0, otherwise the numbers are bound to 'blow up'.

[edit] On second thought, just preset the entire f array to zero (including the boundary elements at x=0 and x=nx), then plug in your initial concentration starting condition on the t=0 row.

 

[nvn]

haitrungdo82: Because you only use the immediately-preceding row of matrix f, then only the immediately-preceding row must preexist, as follows. I numbered the lines, for easy reference.

(1) nt := 2000
(2) nx := 4000
(3) t := 0 ... nt - 1
(4) x := 1 ... nx - 1
(5) f_0,nx := 0
(6) f_0,1 := 1

(7) f_t+1,x := f_t,x - (0.5*v*deltat/deltax)*(f_t,x+1 - f_t,x-1) + (D*deltat/deltax^2)*(f_t,x-1 - 2*f_t,x + f_t,x+1)

Notice, line 5 is the only line you were missing, the absence of which caused the j = 4000 column in matrix f to not exist. Line 6 sets the initial condition, presumably according to your professor (?).

Also, I am not familiar with the formula in line 7, but perhaps you are missing a deltax in the denominator of the second term, which I inserted in line 7? If line 7 crashes, as gneill mentioned, start out with nt and/or nx small, and see how large they can go before line 7 crashes. Then try to figure out what causes it to crash.

 

[berkeman]

Hi qneill,

Btw, if you would like to email me, feel free to do so. My email is << e-mail address deleted by Moderators >>Thank you!

We discourage the posting of personal e-mail addresses in the open forums -- it attracts spam bots. Please use the PM system to initiate private conversations.

 

[haitrungdo82]

Suppose you define your index variables as:

nt := 2000
nx := 4000

t := 0 .. nt - 1
x = 1 .. nx - 1

Your difference equation references f_t+1,x on the LHS, while on the RHS it needs values for f_t,x , f_t,x-1 , f_t,x+1 . This means that it is looking for values calculated at the previous timestep, t = 0, when it's assigning values to f_t+1,__. That's fine.

Now check the x indexes. On the LHS they range from 1 though 3999. On the right they want to go from 0 (for f_t,x-1) through 4000 (for f_t,x+1) at the ends of the range. That means you're going to need to pre-assign values for x=0 and x=4000 for all of t.

Presumably, since you mention that this is a diffusion problem, these pre-assigned values will represent the boundary values for the diffusion field.

I would suggest, then, setting up the entire first row of your f_0,x with the initial concentration pattern, and then setting:

Cb = {some value representing the concentration at the boundaries (could be zero)}
f_t,0 = Cb
f_t,nx = Cb

and then proceed with the diffusion calculations.

You're going to have to work out the issues with the weighting term v*Dt/2. I don't think it should have units or be greater than 1.0, otherwise the numbers are bound to 'blow up'.

[edit] On second thought, just preset the entire f array to zero (including the boundary elements at x=0 and x=nx), then plug in your initial concentration starting condition on the t=0 row.

Hi gnell,

I've just come back from a few day break! Thank you very much for your advice! I think the problem has been solved. Thank you again! So, is it possible to graph this function in 3D, with 3 axis: x, t and f_t,x? I tried using 3D scatter plot. After I put "f" into the placeholder, the program RAN WITHOUT ERRORS, but I saw nothing. I tried changing the axis limits but I still saw nothing. Do you have any suggestion? Thanks!

 

[haitrungdo82]

haitrungdo82: Because you only use the immediately-preceding row of matrix f, then only the immediately-preceding row must preexist, as follows. I numbered the lines, for easy reference.

(1) nt := 2000
(2) nx := 4000
(3) t := 0 ... nt - 1
(4) x := 1 ... nx - 1
(5) f_0,nx := 0
(6) f_0,1 := 1

(7) f_t+1,x := f_t,x - (0.5*v*deltat/deltax)*(f_t,x+1 - f_t,x-1) + (D*deltat/deltax^2)*(f_t,x-1 - 2*f_t,x + f_t,x+1)

Notice, line 5 is the only line you were missing, the absence of which caused the j = 4000 column in matrix f to not exist. Line 6 sets the initial condition, presumably according to your professor (?).

Also, I am not familiar with the formula in line 7, but perhaps you are missing a deltax in the denominator of the second term, which I inserted in line 7? If line 7 crashes, as gneill mentioned, start out with nt and/or nx small, and see how large they can go before line 7 crashes. Then try to figure out what causes it to crash.

Hi nvn,

I have just come back from a short break! You and gnell gave me very useful advice and I think my problem has been solved. You probably have read the question I asked gnell about the graph. Basically, I'm just wondering if it's possible to graph this function in 3D with 3 axis: x, t and f_t,x. I tried using 3D scatter plot and changed the axis limit. Although the program ran and did not show any error, I did not see any line. Could you give me some suggestions? Thank you very much!

 

[gneill]

Hi gnell,

I've just come back from a few day break! Thank you very much for your advice! I think the problem has been solved. Thank you again! So, is it possible to graph this function in 3D, with 3 axis: x, t and f_t,x? I tried using 3D scatter plot. After I put "f" into the placeholder, the program RAN WITHOUT ERRORS, but I saw nothing. I tried changing the axis limits but I still saw nothing. Do you have any suggestion? Thanks!

You should now have array f filled with the data values. Try changing the plot type to "surface plot".

 

[haitrungdo82]

You should now have array f filled with the data values. Try changing the plot type to "surface plot".

Hi gnell,

Interestingly, when I converted the file from MathCAD to pdf, the graph actually appeared. May be the MathCAD program in my computer is not good enough to show the graph. They are attached below. I think it doesn't matter which 3D type I chose, surface or 3D scatter, the graphs looked mostly the same. I'm afraid that they look too dense and doesn't look like the graph from the question. I think may be because delta t is too small, do you think so?

I also attached the graph from the question. I'm not sure if it's possible to make a graph look like this because I don't know whether this graph is just an illustration or an actual graph from MathCAD. Do you have any suggestion? Thank you very much!

 

[haitrungdo82]

You should now have array f filled with the data values. Try changing the plot type to "surface plot".

Hi gnell,

Actually delta t := 0.01 min and delta x := 5m are the best pair of variables that I was able to manage before the program crashed and after including your and nvn's suggestions.

Minh

 

[gneill]

The “Graph from Question” appears to have been created with a drawing program; I doubt that it represents real output from a simulation, but rather indicates the features that are to be expected from a simulation.

Your graphs appear to be so dense because they are! You’re generating a lot of data points from all those iterations. You can probably thin them out for plotting purposes by either building a new array by picking out every nth row and column (for some choice of n), or by just doing some runs with fewer time steps and larger distance steps.

Speaking of step sizes, I think you’ll have to spend a bit of time thinking about what the simulation is supposed to represent, and what the multiplying factors (D and V/2) are supposed to do. Usually in a simulation like this you’ve got a “playing field” of some given dimensions and a desired elapsed time. Both of these are then broken down into time and distance steps that cover the whole field over the duration of the simulation. In this case you’ve got a 1D space (a line) where some values propagate along it time-step by time-step. Without some idea of how the simulation is supposed to map onto the real physics, it will be difficult to know what parameter sizes are reasonable.

So, how long is this line supposed to be in “real life”? Is it a line a meter long? A kilometer? How about 100 kilometers? I can guess that the latter is more likely given that the suggested velocity parameter is given in kilometers per hour. The D parameter seems to be specified in m² per minute, suggesting that a suitable time span for the simulation would be on the order of (simulated) hours or days.

Once you’ve established the dimensions of the simulation (time and space), and the numbers of steps of each you want to take (granularity), then that sets the step sizes. For example, say that the line represents 100km and the time is to span 12 hours. If nx is the number of special steps and nt the number of time steps, the step sizes are $$\Delta$$t = 12hr/nt and $$\Delta$$x = 100km/nx.

There’s not much more I can help with regarding the simulation, since I don’t know the context or have knowledge of the course material.

 

[haitrungdo82]

The “Graph from Question” appears to have been created with a drawing program; I doubt that it represents real output from a simulation, but rather indicates the features that are to be expected from a simulation.

Your graphs appear to be so dense because they are! You’re generating a lot of data points from all those iterations. You can probably thin them out for plotting purposes by either building a new array by picking out every nth row and column (for some choice of n), or by just doing some runs with fewer time steps and larger distance steps.

Speaking of step sizes, I think you’ll have to spend a bit of time thinking about what the simulation is supposed to represent, and what the multiplying factors (D and V/2) are supposed to do. Usually in a simulation like this you’ve got a “playing field” of some given dimensions and a desired elapsed time. Both of these are then broken down into time and distance steps that cover the whole field over the duration of the simulation. In this case you’ve got a 1D space (a line) where some values propagate along it time-step by time-step. Without some idea of how the simulation is supposed to map onto the real physics, it will be difficult to know what parameter sizes are reasonable.

So, how long is this line supposed to be in “real life”? Is it a line a meter long? A kilometer? How about 100 kilometers? I can guess that the latter is more likely given that the suggested velocity parameter is given in kilometers per hour. The D parameter seems to be specified in m² per minute, suggesting that a suitable time span for the simulation would be on the order of (simulated) hours or days.

Once you’ve established the dimensions of the simulation (time and space), and the numbers of steps of each you want to take (granularity), then that sets the step sizes. For example, say that the line represents 100km and the time is to span 12 hours. If nx is the number of special steps and nt the number of time steps, the step sizes are $$\Delta$$t = 12hr/nt and $$\Delta$$x = 100km/nx.

There’s not much more I can help with regarding the simulation, since I don’t know the context or have knowledge of the course material.

Hi gnell,

Can the probability density be ever greater than 1? I think it can not exceed 1. But as the way the program is written, it can be very very large! I'm confused!

 

[gneill]

Hi gnell,

Can the probability density be ever greater than 1? I think it can not exceed 1. But as the way the program is written, it can be very very large! I'm confused!

Usually, when things go "unphysical" in a simulation it's an indication that something has gone wrong with the basic assumptions, such as unrealistic parameters for the given model, or steps too big so that they step over and miss important interactions, or too small so that error accumulates and rivals the "signal".

I think you'll have to do some investigation into what it really is that is being modeled so that all the parameters can be chosen 'realistically'.

 

[haitrungdo82]

Usually, when things go "unphysical" in a simulation it's an indication that something has gone wrong with the basic assumptions, such as unrealistic parameters for the given model, or steps too big so that they step over and miss important interactions, or too small so that error accumulates and rivals the "signal".

I think you'll have to do some investigation into what it really is that is being modeled so that all the parameters can be chosen 'realistically'.

Hi gnell,

If you feel that's enough, you can ignore my question! I think the simulation is in shape now except some minor strange occurences (I believe!).

If we disregard the graph (which is too dense to be visible) and just look at the table of values, the data are actually following an expected trend in which the smoke stack dispenses toward left and right sides from the initial point (f(t=0, x=2000)) as time goes on, just like the illustration image.

As you look at the attached MathCAD sheet below, all the values are now less than 1. But there are negative values?? From my observation, the negative and positive values appear in an alternating fashion. Is it a typical problem encountered in numerial simulation?

I'm trying to vary the step sizes and the number of steps. But I'm just wondering if you have any better suggestions.

Thanks a lot, gnell!

 

[gneill]

You'll have to decide if negative values have proper physical meaning for what is being modeled. If not, there is some sort of problem, either with the model itself, the way it is set up (parameters), or the interpretation of the results.

If this is a well-known model, then there are probably papers written which discuss its use, stability, and overall behavior.

 

[nvn]

haitrungdo82: No matter what you do, the concentration of the pollutant rapidly decays. Only a single glob of your pollutant is being released (at t = 0), from a point source (x = 1); therefore, f_0,1 := 1. Just the slightest wind will cause this pollutant to travel in the positive x direction. And the diffusion will cause the pollutant to slightly spread out in the positive and negative x directions, but not as fast as the wind pushes the pollutant in the positive x direction. Therefore, your problem is one-sided, meaning all pollutant travels from x = 1 toward the positive x axis. I.e., virtually no pollutant travels in the negative x direction. Therefore, you can start your problem at x = 1, meaning you can let the pollutant release occur at x = 1. Thus, f_0,1 := 1.

I recommend you try the following.

(1) deltat := 0.035 s
(2) deltax := 0.0002 m
(3) v := 0.003 m/s
(4) D = 0.000 017 (m^2)/min
(5) nt := 90
(6) nx := 60
(7) t := 0 ... nt - 1
(8) x := 1 ... nx - 1
(9) f_0,nx := 0
(10) f_0,1 := 1
(11) (same as line 7 in post 9)

But when you plot the 3-D scatter plot for this data set, I recommend you do not plot the first ten rows of matrix f. This will improve your visibility of the plot.

Note that if you run the problem out too far, or use parameter values too large or small, it would produce negative values, and garbage, in matrix f. In reality, perhaps you would need to put a conditional statement that causes f < 0 to be changed to f = 0; but that would use up or waste execution time. And I doubt it would eliminate all garbage. The minimum possible value of the pollutant concentration (or probability of concentration) is zero, not negative. And the maximum value is f = 1.

Now change v to 0.0005 m/s, and you will see that the pollutant does not move along the positive x-axis as rapidly, during this ~3 s simulation.

Only a single puff of your pollutant is released, which occurs at t = 0. If you changed your initial condition to f_t,0 := 1, f_0,0 := 0, then your pollutant would be continuously emitted from the pollution source. Then you would see just a wall of smoke moving along the x axis, which is somewhat less eventful.

Please let us know if you uncover anything else about the diffusion equation, or how to use it. I am currently just groping through the dark.

 

[haitrungdo82]

haitrungdo82: No matter what you do, the concentration of the pollutant rapidly decays. Only a single glob of your pollutant is being released (at t = 0), from a point source (x = 1); therefore, f_0,1 := 1. Just the slightest wind will cause this pollutant to travel in the positive x direction. And the diffusion will cause the pollutant to slightly spread out in the positive and negative x directions, but not as fast as the wind pushes the pollutant in the positive x direction. Therefore, your problem is one-sided, meaning all pollutant travels from x = 1 toward the positive x axis. I.e., virtually no pollutant travels in the negative x direction. Therefore, you can start your problem at x = 1, meaning you can let the pollutant release occur at x = 1. Thus, f_0,1 := 1.

I recommend you try the following.

(1) deltat := 0.035 s
(2) deltax := 0.0002 m
(3) v := 0.003 m/s
(4) D = 0.000 017 (m^2)/min
(5) nt := 90
(6) nx := 60
(7) t := 0 ... nt - 1
(8) x := 1 ... nx - 1
(9) f_0,nx := 0
(10) f_0,1 := 1
(11) (same as line 7 in post 9)

But when you plot the 3-D scatter plot for this data set, I recommend you do not plot the first ten rows of matrix f. This will improve your visibility of the plot.

Note that if you run the problem out too far, or use parameter values too large or small, it would produce negative values, and garbage, in matrix f. In reality, perhaps you would need to put a conditional statement that causes f < 0 to be changed to f = 0; but that would use up or waste execution time. And I doubt it would eliminate all garbage. The minimum possible value of the pollutant concentration (or probability of concentration) is zero, not negative. And the maximum value is f = 1.

Now change v to 0.0005 m/s, and you will see that the pollutant does not move along the positive x-axis as rapidly, during this ~3 s simulation.

Only a single puff of your pollutant is released, which occurs at t = 0. If you changed your initial condition to f_t,0 := 1, f_0,0 := 0, then your pollutant would be continuously emitted from the pollution source. Then you would see just a wall of smoke moving along the x axis, which is somewhat less eventful.

Please let us know if you uncover anything else about the diffusion equation, or how to use it. I am currently just groping through the dark.

Hi nvn,

Right! Below is a direct quote from the question:

"A stack emits at time t = 0 a short burst of pulluting gas that is carried by a steady wind (v = vx = 10 mi/h) (in x direction) while it undergoes lateral dispersion (in x direction). Stochastic dispersion is a process similar to diffusion but under the influence of a driving field..."

I forgot the fact that the wind actually blew only in one direction. So does it mean the illustration is misleading? If you look at one of the above attachments, named "Graph from Question", which I copied directly from the question, it shows a Gaussian distribution at every time step t. I'm not sure if my Professor included this figure intentionally or by mistake.

In one of his hints, he said that if we ran into numerical problems manifesting themselves by unphysical fluctuations of the calculated probability distribution, increase the number of bins in t and x, which increased the accuracy. May be he mentioned the problem that I am encountering right now. I did it, but MathCAD had a limit! When I increased t or x up to 6000 or 7000, MathCAD’s memory was not able to perform the computation.

But I’ll take your advice because, right now, my values are either positive or negative and several probability densities are even greater than 1, which are very unphysical

Thank you very much for your advice!

 

[haitrungdo82]

haitrungdo82: No matter what you do, the concentration of the pollutant rapidly decays. Only a single glob of your pollutant is being released (at t = 0), from a point source (x = 1); therefore, f_0,1 := 1. Just the slightest wind will cause this pollutant to travel in the positive x direction. And the diffusion will cause the pollutant to slightly spread out in the positive and negative x directions, but not as fast as the wind pushes the pollutant in the positive x direction. Therefore, your problem is one-sided, meaning all pollutant travels from x = 1 toward the positive x axis. I.e., virtually no pollutant travels in the negative x direction. Therefore, you can start your problem at x = 1, meaning you can let the pollutant release occur at x = 1. Thus, f_0,1 := 1.

I recommend you try the following.

(1) deltat := 0.035 s
(2) deltax := 0.0002 m
(3) v := 0.003 m/s
(4) D = 0.000 017 (m^2)/min
(5) nt := 90
(6) nx := 60
(7) t := 0 ... nt - 1
(8) x := 1 ... nx - 1
(9) f_0,nx := 0
(10) f_0,1 := 1
(11) (same as line 7 in post 9)

But when you plot the 3-D scatter plot for this data set, I recommend you do not plot the first ten rows of matrix f. This will improve your visibility of the plot.

Note that if you run the problem out too far, or use parameter values too large or small, it would produce negative values, and garbage, in matrix f. In reality, perhaps you would need to put a conditional statement that causes f < 0 to be changed to f = 0; but that would use up or waste execution time. And I doubt it would eliminate all garbage. The minimum possible value of the pollutant concentration (or probability of concentration) is zero, not negative. And the maximum value is f = 1.

Now change v to 0.0005 m/s, and you will see that the pollutant does not move along the positive x-axis as rapidly, during this ~3 s simulation.

Only a single puff of your pollutant is released, which occurs at t = 0. If you changed your initial condition to f_t,0 := 1, f_0,0 := 0, then your pollutant would be continuously emitted from the pollution source. Then you would see just a wall of smoke moving along the x axis, which is somewhat less eventful.

Please let us know if you uncover anything else about the diffusion equation, or how to use it. I am currently just groping through the dark.

Hi again, nvn,

I just forgot! The question actually required me to inspect the evolution of the dispersion of the pollutant over a lateral distance -10km <= x-x0 <= 10km and a time span of 200 min. That's why I had to include a large number of t and x. And if the negative x direction can be ignored, why would he mind asking me to inspect the dispersion up to -10 km?

I also used "if" to try to eliminate all negative values. It worked! But...another problem occured. My numerical model went from Gaussian distribution to a "disorder" distribution. So, I was afraid that if the negeative values actually affect the overall behavior, so changing them to positive may incorectly alter the model. But...the values should not be negative...So,...confusing!

 

[nvn]

haitrungdo82: Where did the diffusion coefficient value, D = 0.000 017 (m^2)/min, come from? Is that given in the question?

 

[haitrungdo82]

haitrungdo82: Where did the diffusion coefficient value, D = 0.000 017 (m^2)/min, come from? Is that given in the question?

Hi nvn,

Yes. All original data were given in the question. But, t, x, delta t, and delta x can be varied. However, I think, eventually, the program should be large enough so that I am able to inspect the model behavior in the range [-10km,10km] for 200 mins.

I'm still thinking about your suggestion. It makes sense to me. But it seems like my Professor didn't think so. Why did he want to see the behavior as far as -10km if the values can be ignored. Do you agree?

 

[haitrungdo82]

haitrungdo82: Where did the diffusion coefficient value, D = 0.000 017 (m^2)/min, come from? Is that given in the question?

Hi nvn,

Also, since -10km <= x-x0 <= 10 km, probably x0 is the center of the range, which is probably x0 = 2000 if nx = 4000

 

[nvn]

haitrungdo82: We are still assimilating the information you posted today, so give us quite some time to think it over. But first, even before we have uncovered anything yet, I think we can say the current diffusion coefficient must be grossly incorrect. Shouldn't it be at least something like D = 0.000 017 (km^2)/min, or more, instead of 0.000 017 (m^2)/min? Could someone let us know?

haitrungdo82 (or anyone), do you have a table of typical diffusion coefficients for something like typical smoke in air? See if you can research this, to see if there is a typographic mistake in the current D value or units.

haitrungdo82, do you have any reason to believe any parameter value, or units, in the given question could be a typographic mistake? Or are you certain the professor did not make any mistake?

Don't worry about the -10 km just yet. We can address that later. We can plot to any negative value.

 

[haitrungdo82]

haitrungdo82: We are still assimilating the information you posted today, so give us quite some time to think it over. But first, even before we have uncovered anything yet, I think we can say the current diffusion coefficient must be grossly incorrect. Shouldn't it be at least something like D = 0.000 017 (km^2)/min, or more, instead of 0.000 017 (m^2)/min? Could someone let us know?

haitrungdo82 (or anyone), do you have a table of typical diffusion coefficients for something like typical smoke in air? See if you can research this, to see if there is a typographic mistake in the current D value or units.

haitrungdo82, do you have any reason to believe any parameter value, or units, in the given question could be a typographic mistake? Or are you certain the professor did not make any mistake?

Don't worry about the -10 km just yet. We can address that later. We can plot to any negative value.

Hi nvn,

Thank you very much for your kindness! This diffusion constant is actually a corrected version. He might have made a mistake again, but it should be rare because he probably ran the program and figured out the mistake with the original version, then corrected it, and the program should have run successfully with this corrected coefficient.

I'm just wondering how important the difussion constant is in this case. Let's say this is an incorrect coefficient. Can we still construct the model without thinking about it? How is it related to my problem, specifically, the unphysical negative probabilities? I actually tried changing the magnitude of the coefficient, but it didn't seem to change the overall behavior of the model. Am I correct?

OK, if the coefficient is really important, can you approximate a reasonable value? As long as it makes sense, we can change everything and comment on it. It is also part of this assignment

nvn, keep me updated with your new ideas. I probably have to go on with this model and submit it, no matter how bad it is, on the Monday of next week.

But, again, keep me up to date, even at the last minute. This model is hard to construct but easy to change. Thank you very much!

 

[nvn]

haitrungdo82: I agree with all of your statements in all of your posts yesterday, except for the 200 min. An analysis time of 30 min is better. Your wind velocity is 4.47 m/s. If you run the simulation more than, say, 30 min, the pollutant has already been moved off of the 10 km analysis space by the wind. Therefore, it is irrelevant to run the simulation longer.

D = 0.000 017 (m^2)/min appears incorrect. A value this small would have no effect on the problem, and would be irrelevant.

First, just for fun, try the following.

(1) nt := 2000
(2) nx := 200
(3) deltat := 6 s
(4) deltax := 50 m
(5) v := 0 m/s
(6) D := 0.01 (km^2)/min
(7) f_0,nx := 0
(8) f_0,0.5*nx := 1

Then try the following.

(9) nt := 2000
(10) nx := 2000
(11) deltat := 0.90 s
(12) deltax := 10 m
(13) v := 4.470 m/s
(14) D := 0.0017 (km^2)/min
(15) f_0,nx := 0
(16) f_0,0.5*nx := 1

You could also try the following, but it does not appear to be better than lines 9 through 16. Therefore, lines 9 through 16 seem better.

(17) nt := 3600
(18) nx := 2200
(19) deltat := 0.50 s
(20) deltax := 9 m
(21) v := 4.470 m/s
(22) D := 0.0017 (km^2)/min
(23) f_0,nx := 0
(24) f_0,0.5*nx := 1

 

[haitrungdo82]

haitrungdo82: I agree with all of your statements in all of your posts yesterday, except for the 200 min. An analysis time of 30 min is better. Your wind velocity is 4.47 m/s. If you run the simulation more than, say, 30 min, the pollutant has already been moved off of the 10 km analysis space by the wind. Therefore, it is irrelevant to run the simulation longer.

D = 0.000 017 (m^2)/min appears incorrect. A value this small would have no effect on the problem, and would be irrelevant.

First, just for fun, try the following.

(1) nt := 2000
(2) nx := 200
(3) deltat := 6 s
(4) deltax := 50 m
(5) v := 0 m/s
(6) D := 0.01 (km^2)/min
(7) f_0,nx := 0
(8) f_0,0.5*nx := 1

Then try the following.

(9) nt := 2000
(10) nx := 2000
(11) deltat := 0.90 s
(12) deltax := 10 m
(13) v := 4.470 m/s
(14) D := 0.0017 (km^2)/min
(15) f_0,nx := 0
(16) f_0,0.5*nx := 1

You could also try the following, but it does not appear to be better than lines 9 through 16. Therefore, lines 9 through 16 seem better.

(17) nt := 3600
(18) nx := 2200
(19) deltat := 0.50 s
(20) deltax := 9 m
(21) v := 4.470 m/s
(22) D := 0.0017 (km^2)/min
(23) f_0,nx := 0
(24) f_0,0.5*nx := 1

Hi nvn,

I just read your post. I think I got what you mean! Let's say we are faithful with my Professor's assumption, t= 200 min. Then, we take a "screen shot" from -10km to 10km, and ignore the rest. I guess you meant after 200 min, most of the pollutants already traveled far away from 10 km and the probability is almost 0 throughout the area of our inspection. Am I right?

It sounds reasonable! I'll take your advice and let you know what I think. Be in touch!

Thank you!

 

[haitrungdo82]

haitrungdo82: I agree with all of your statements in all of your posts yesterday, except for the 200 min. An analysis time of 30 min is better. Your wind velocity is 4.47 m/s. If you run the simulation more than, say, 30 min, the pollutant has already been moved off of the 10 km analysis space by the wind. Therefore, it is irrelevant to run the simulation longer.

D = 0.000 017 (m^2)/min appears incorrect. A value this small would have no effect on the problem, and would be irrelevant.

First, just for fun, try the following.

(1) nt := 2000
(2) nx := 200
(3) deltat := 6 s
(4) deltax := 50 m
(5) v := 0 m/s
(6) D := 0.01 (km^2)/min
(7) f_0,nx := 0
(8) f_0,0.5*nx := 1

Then try the following.

(9) nt := 2000
(10) nx := 2000
(11) deltat := 0.90 s
(12) deltax := 10 m
(13) v := 4.470 m/s
(14) D := 0.0017 (km^2)/min
(15) f_0,nx := 0
(16) f_0,0.5*nx := 1

You could also try the following, but it does not appear to be better than lines 9 through 16. Therefore, lines 9 through 16 seem better.

(17) nt := 3600
(18) nx := 2200
(19) deltat := 0.50 s
(20) deltax := 9 m
(21) v := 4.470 m/s
(22) D := 0.0017 (km^2)/min
(23) f_0,nx := 0
(24) f_0,0.5*nx := 1

Hi nvn,

Thumb up! All the negative values disappeared and the model is physically perfect! It seems like the diffusion coefficient is the cause of all messy, unphysical things. I'm just curious! What were you thinking when you were trying to eliminate the negative values? I'm very curious about this! Thank you!

 

[nvn]

haitrungdo82: My thinking was very vague. We know 0.5*(f_t,x+1 - f_t,x-1)/deltax is the first-order finite difference approximation of the first derivative, if I recall correctly. And v*deltat is an x (horizontal) distance. Therefore, this horizontal distance times the above-mentioned first derivative (which would be the tangent of the slope angle) should be a vertical (y) distance, if I recall correctly.

We know (f_t,x-1 - 2*f_t,x + f_t,x+1)/deltax^2 is the first-order finite central difference approximation of the second derivative (curvature), if I recall correctly. Therefore, similarly, when we multiply this value times D*deltat, we get another y distance. I do not recall the exact details.

Therefore, I vaguely envisioned that if we supply a D value too large or too small, the equation was overshooting below y = 0, or overshooting in the horizontal direction (?), which, in either case, would give garbage. I envisioned that if D is too small, it has no effect on the dispersion. If D is too large, the equation overshoots. I envisioned if deltax is too small, we get almost a singularity, called overflow. We want to increase D, but I envisioned that we cannot increase it too much unless we decrease deltat; otherwise, the problem will overshoot, as explained above. I envisioned that you cannot decrease deltat below a certain value, or else the problem will crash or overshoot. And I envisioned that you cannot decrease deltax below a certain value, or else the problem will crash or overshoot. E.g., if we want to simulate out to 10 km, it seems deltat cannot be less than 0.5 s, and deltax cannot be less than ~8 m. I do not understand very well what was going on, nor how to control the diffusion equation. If you uncover anything else, feel free to let us know.

I had to increase D seemingly too much (?); but that is the only way I could get the problem to run. In your given problem, D is overpowered by the relatively high wind velocity.

In the "Graph from Question" file, we see a lot of diffusion, and almost no wind. We can see that the centerline of the probability density moves slightly in the negative x direction, which means there is a very, very slight wind in the negative x direction.

By the way, to reply to a post, did you know you can just press the "New Reply" button? You don't need to press the "Quote" button on every reply. Also, always leave a space between a numeric value and its following unit symbol. E.g., 10 km, not 10km. See the international standard for writing units (ISO 31-0).

If you uncover anything else, feel free to let us know.

 

[haitrungdo82]

haitrungdo82: My thinking was very vague. We know 0.5*(f_t,x+1 - f_t,x-1)/deltax is the first-order finite difference approximation of the first derivative, if I recall correctly. And v*deltat is an x (horizontal) distance. Therefore, this horizontal distance times the above-mentioned first derivative (which would be the tangent of the slope angle) should be a vertical (y) distance, if I recall correctly.

We know (f_t,x-1 - 2*f_t,x + f_t,x+1)/deltax^2 is the first-order finite central difference approximation of the second derivative (curvature), if I recall correctly. Therefore, similarly, when we multiply this value times D*deltat, we get another y distance. I do not recall the exact details.

Therefore, I vaguely envisioned that if we supply a D value too large or too small, the equation was overshooting below y = 0, or overshooting in the horizontal direction (?), which, in either case, would give garbage. I envisioned that if D is too small, it has no effect on the dispersion. If D is too large, the equation overshoots. I envisioned if deltax is too small, we get almost a singularity, called overflow. We want to increase D, but I envisioned that we cannot increase it too much unless we decrease deltat; otherwise, the problem will overshoot, as explained above. I envisioned that you cannot decrease deltat below a certain value, or else the problem will crash or overshoot. And I envisioned that you cannot decrease deltax below a certain value, or else the problem will crash or overshoot. E.g., if we want to simulate out to 10 km, it seems deltat cannot be less than 0.5 s, and deltax cannot be less than ~8 m. I do not understand very well what was going on, nor how to control the diffusion equation. If you uncover anything else, feel free to let us know.

I had to increase D seemingly too much (?); but that is the only way I could get the problem to run. In your given problem, D is overpowered by the relatively high wind velocity.

In the "Graph from Question" file, we see a lot of diffusion, and almost no wind. We can see that the centerline of the probability density moves slightly in the negative x direction, which means there is a very, very slight wind in the negative x direction.

By the way, to reply to a post, did you know you can just press the "New Reply" button? You don't need to press the "Quote" button on every reply. Also, always leave a space between a numeric value and its following unit symbol. E.g., 10 km, not 10km. See the international standard for writing units (ISO 31-0).

If you uncover anything else, feel free to let us know.

Hi nvn,

Thank you very much for your advice! I was afraid that if I did not press the "Quote" button, you would not be aware of my new messages. If this is not the case, I will press the "New Reply" next time.

This is probably my last question for you. I was going to take the natural log (ln) of each value in the array. However, since there are many "0"s, the program will crash. There are 2 ways I can solve this problem. I can either convert all "0" into "1" or skip the "0" values. I don't know how to skip some selected values in the array, so I chose the first method.

I used "if" to convert all "0" to "1" as you can see from the attachment. However, only the upper part of the array was actually converted. The rest still contained "0". I copied 2 different areas of the same array in the attachment to illustrate my point.

So, my question is whether I missed anything when I applied the condition. Or if you know how to skip some values in an array, it would be great and I would like to learn.

Thank you very much for your help

 

[haitrungdo82]

You'll have to decide if negative values have proper physical meaning for what is being modeled. If not, there is some sort of problem, either with the model itself, the way it is set up (parameters), or the interpretation of the results.

If this is a well-known model, then there are probably papers written which discuss its use, stability, and overall behavior.

Hi gnell,

I just posted a question regarding how to convert some values in an array to other values or how to skip some values in an array. Hope you or nvn can help me since both of you are very good at MathCAD. Thank you very much!

 

[nvn]

haitrungdo82: We will indeed be aware of your new messages if you press "New Reply."

Try changing your program to something like f5_t,x := {f5_t,x if f5_t,x > 10^-8; 1 otherwise}, and see if it then works. You can adjust the -8 exponent to any value you prefer, such as -6, -8, -10, etc. Also, ensure you are displaying enough decimal places, or scientific notation, in your matrix to enable you to see the actual values stored in your matrix.