复件 frozen_wave.cpp

数值模拟交通流模型Payne-weber压缩流体

/*	Aaron Bernard
Math 164: Scientific Computing

This program uses the Frozen-Wave Model as a basis for numerical computation
of a hardcoded PDE.  It does not support input/output.  Matrices are
hardcoded through #define statements.  Two good features to fix in this code
would be input of parameters, as well as use of STL to get rid of hardcoding
in old C style matrices.

The numerical method used to compute the PDE is the upwind method, or
Godunov method.  This is described in a paper by Jin and Zhang called
"Studies on a Nonequilibrium Continuum Traffic Flow Model with Frozen
'Sound Wave' Speed" (2002).  A hyperlink to the site is:

http://www.its.berkeley.edu/publications/ITSReviewonline/spring2003/trb2003/jin.pdf
*/

#include <iostream>
#include <fstream>
#include <cmath>
#include <assert.h>
using namespace std;

#define xsteps 100
#define tsteps 100
#define tint 10
#define xint 10

double mid(double left, double right);

int main(int argc, char* argv[])
{
	// The following two components are part of the vector u
	double p[xsteps][tsteps] = {0};
	double v[xsteps][tsteps] = {0};

	double dx = double(xint)/xsteps;
	double dt = double(tint)/tsteps;
	double isl = dt/dx;			// Inverse slope
	double e = 0.01;			// epsilon
	double c = 0.5;				// mean free speed
	double m = 0.0;				// used for midpoint

	fstream fout;
	fout.open("upwind.m",ios_base::out);

	// Define initial conditions for u
	for (int i = 0; i < xsteps; ++i) {
		p[i][0] = (1.0/xsteps)*i + e;
		v[i][0] = 1 - (1.0/xsteps)*i + e;
		fout << p[i][0] << " " << v[i][0] << endl;
	}

	// Evolve the solution in time according to the index j
	for (int j = 1; j < tsteps - 1; ++j) {
		// Define boundary condition of x = 0 off of previous timestep
		m = mid(v[0][j-1],v[1][j-1]) - c;
		p[0][j] = (m > 0 ? p[0][j-1] : mid(p[0][j-1],p[1][j-1]));
		v[0][j] = (m > 0 ? v[0][j-1] : m);

		// Calculate the solution of u = [p,v] in this, the j+1 th timestep
		for (int i = 1; i < xsteps - 1; ++i) {		
			p[i][j+1] = p[i][j] + isl*(	// Original Point + Slope Times Sum of:
				mid(p[i][j-1],p[i-1][j-1])*mid(v[i][j-1],v[i-1][j-1]) -	// Lo Flux
				mid(p[i][j-1],p[i+1][j-1])*mid(v[i][j-1],v[i+1][j-1])   // Hi Flux
				);

			v[i][j+1] = v[i][j] + isl*(	// Original Point + Slope Times Lo,Hi Flux
				(.5*pow(mid(v[i][j-1],v[i-1][j-1]),2) - c*mid(v[i][j-1],v[i-1][j-1])) - 
				(.5*pow(mid(v[i][j-1],v[i+1][j-1]),2) - c*mid(v[i][j-1],v[i+1][j-1]))
				);
		}} // End of both for loops


	// Output the solution data
	fout << "[";
	for (int j = 0; j < tsteps; ++j) {
		for (int i = 0; i < xsteps; ++i) {
			fout << " " << p[i][j];
		}
		fout << endl;
	}
	fout << "]";
	fout.close();

	return 0;
}

/* This function is used to create an average value between two cells
e.g. U^i_(j+i/2) = (U^i_j + U^i_j+1)/2 */
double mid(double left, double right)
{
	double v = (left+right)/2.0;
	return v;
}