Appendix B:
Source Code for fwdeuler.cpp

//Implementation of explicit Euler's method
//CDT Zac Miller, 11-20 April 1999
#include <ctype.h>
#include <iostream.h>
#include <math.h>
#include "vector.h"
#define LD double
//constants
const int max_points = 51;
const LD film_level = -.5;
 
//function definitions
LD S1 (LD);
LD S2 (LD);
LD drhodx(LD, vector, vector, vector);
vector dudx(vector[], int, LD);
void file_input (int&, LD&, LD&, LD [], vector [], vector [], vector []);
void file_output (vector[], int);
//S1
LD S1 (LD x)
{
return (20 * x - 10);
}
//S2
LD S2 (LD x)
{
return (4 * pow((x - .5),2));
}
 
 
 
//--------------------------------------
//| drho/dx |
//--------------------------------------
LD drhodx(LD rho, vector S, vector U, vector dudx)
{
LD result;
result = rho * dot(S, dudx) /
(normal(subtract(S, multiply(rho, U))) - dot(S, U) + rho);
return (result);
}
 
//--------------------------------------
//| du/dx (linear approximation) |
//--------------------------------------
vector dudx(vector u[], int i, LD deltax)
{
vector result;
if (i == 0)
result = multiply(1/deltax, subtract(u[i+1], u[i]));
else
result = multiply(1/deltax, subtract(u[i], u[i-1]));
return (result);
}
 
//--------------------------------------
//| Main program |
//--------------------------------------
main()
{
int num_points = 0, rho_initial_x_i;
LD x[max_points], rho[max_points], deltax, initial_condition;
vector S[max_points], F[max_points], P[max_points], U[max_points];
char flag, junk[10];
file_input (num_points, deltax, initial_condition, x, F, S, U);
rho_initial_x_i = (num_points-1)/2;
rho[rho_initial_x_i] = initial_condition;
P[rho_initial_x_i].a = 0;
P[rho_initial_x_i].b = initial_condition;
for (int i = rho_initial_x_i + 1; i < num_points; i++)
{
rho[i] = rho[i-1] + deltax *
drhodx(rho[i-1], S[i-1], U[i-1], dudx(U, i, deltax));
P[i] = multiply(rho[i], U[i]);
}
for (i = rho_initial_x_i - 1; i >= 0; i--)
{
rho[i] = rho[i+1] - deltax *
drhodx(rho[i+1], S[i+1], U[i+1], dudx(U, i, deltax));
P[i] = multiply(rho[i], U[i]);
}
cout << "\n\n\n\n\n\nDo you want file output? [Y/N] ";
cin.get(flag);
cin.getline(junk,10);
if ((flag == 'Y') || (flag == 'y'))
file_output(P, num_points);
return (0); //end main
}
 
 
 
//--------------------------------------
//| File input module |
//--------------------------------------
void file_input (int& num_points, LD& deltax, LD &initial_condition,
LD x[], vector F[], vector S[], vector U[])
{
//function declaration
//variables
char filename[max_len], data[max_len]; //filename specified by user
FILE *ptr; // Pointer to the file.
LD startx = 0;
 
// Open the file - no error checking done
cout << "\nEuler's Method (explicit) ";
cout << "\n Please enter the filename for processing: ";
cin.getline(filename, max_len);
ptr = fopen(filename,"r");
fgets(data, max_len, ptr);
num_points = int(extract(data, max_len));
fgets(data, max_len, ptr);
startx = extract(data, max_len);
fgets(data, max_len, ptr);
deltax = extract(data, max_len);
fgets(data, max_len, ptr);
initial_condition = extract(data, max_len);
 
for (int i = 0; (i < num_points) && (feof(ptr) == 0); i++)
{
fgets(data, max_len, ptr); // Read next record
F[i].a = extract(data, max_len);
F[i].b = film_level;
}
fclose(ptr); // Close the file.
//calculate all values
for (i = 0; i < num_points; i++)
{
if (i == 0)
x[i] = startx;
else
x[i] = x[i-1] + deltax;
S[i].a = S1(x[i]);
S[i].b = S2(x[i]);
assign(U[i],multiply(-1/normal(F[i]),F[i]));
}
return;
}
 
 
//------------------------------
//|file output |
//------------------------------
void file_output (vector P[], int num_points)
{
char filename[max_len]; //filename specified by user
FILE *fp;
cout << "Please enter the name of the output file: ";
cin.getline(filename, max_len);
fp = fopen(filename,"w"); //open file
fprintf(fp, "\nTable of P values\n");
fprintf(fp, "\n-----------------------------------------------------");
fprintf(fp, "\n\n i \t\t P1 \t\t P2\n");
for (int i = 0; i < num_points; i++)
{
fprintf(fp, "%2d \t\t%+.10f \t\t%+.10f\n", i, P[i].a, P[i].b);
}
fclose(fp); // Closes the output file
}