Code:
In a total population of n individuals there are, at any time t,
x(t) members of the population who are susceptible to a contagious disease,
y(t) infectious carriers of this contagious disease,
z(t) individuals who are recovered and immune.
It is clear x(t) + y(t) + z(t) = n for all non-negative t
The following systems of equations model the epidemic
dx/dt = - lambda x(t)y(t)
dy/dt = lambda x(t)y(t) - mu y(t)
dz/dt = mu y(t)
The first equation says the rate of change in the susceptible population is
deceasing by the infection rate times the contacts measured by $xy$.
The last equation says the rate of change in the recovered population is
increasing by the recovery rate times the number of infected.
The middle equation, says the infection rate is gaining from susceptibles
and losing to the recovered and it is balanced so that population remains constant.
Input a text file of parameters in the order
lambda mu x(0) y(0) z(0) tinit tfinal tdelta
The differential equation must be solved numerically using the numerical
method described below called the improved Euler method Euler's method for dx/dt = f(t,x)
t_new = t_old + delta_t
dxdt_est = f(t_old, x_old)
x_new = x_old + dxdt_est*delta_t
One can improve the estimate by doing
dxdt_est2 = (dxdt_est + f(t_new, x_new))/2
x_new = x_old + dxdt_est2*delta_t
And you need to repeat this four times.
(Note that x_new changes and dxdt_est2 changes but dxdt_est does not change.
The results this solution must be stored in a 3 x n matrix called solution
and the transpose of this matrix must be output to a file.
Also required is a graph ploting x, y, z and dy/dt on the same graph.
Then for each t, number of sick days the integral of y from 0 to t must be
computed by simpson rule and plotted. Since Simpson's rule requires an even
number of data points, interpolate the data points on either side of those times with an odd.
In most be modular, with a main program that calls functions and/or procedures to modular cores.
Communtication cannot be by global variables, input and output needs to be done by
parameters to functions and/or procedures.
Variable names need to be meaningful.
The Code must be Commented.
The code must be correct, and part of your project needs to state how you checked correctness.
Each project requires the use of the 2D array data type explicitly.
Each project requires reading input data from a text file and outputing a graph (and possibly more).
While generally not an issue, a program can be returned to replace the use of an inefficient algorithm with a better one.
Here is the code I already have (the input file is attached)
Code:
#include <stdio.h>
float dxdt (float lambda, float y_new, float x_new); //function for dx/dt = - lambda x(t)y(t)
float dydt (float lambda, float x_new, float y_new, float mu); //function for dy/dt = lambda x(t)y(t) - mu y(t)
float dzdt (float mu, float y_new); //function for dz/dt = mu y(t)
float simp_y (float t_new, float y_new4); //function for integral of y using Simpson's Rule
int main()
{
float lambda; //infection rate
float mu; //recovery rate
float x0; //initial susceptible population, 900
float y0; //initial infected population, 10
float z0; //initial recovered population, 90
float tinit; //intial time, 0 days
float tfinal; //final time, 15 days
float tdelta; //time steps, 1/10th of a day, 0.1
float x_new; //variables for new x, y, and z values
float y_new; //
float z_new; //
float x_new2; //variables for new x, y, and z values in succesive improvements
float y_new2; //
float z_new2; //
//
float x_new3; //
float y_new3; //
float z_new3; //
//
float x_new4; //
float y_new4; //
float z_new4; //
float t_new; //variable for new time value
float t_old; //variable for old time value
float n; //total number of steps, 10 steps/day * 15 days = 150 steps
float dxdt_est; //1st estimate improvement for dxdt, dydt, and dzdt
float dydt_est; //
float dzdt_est; //
float dxdt_est2; //2nd estimate improvement for dxdt, dydt, and dzdt
float dydt_est2; //
float dzdt_est2; //
float dxdt_est3; //3rd estimate improvement for dxdt, dydt, and dzdt
float dydt_est3; //
float dzdt_est3; //
float dxdt_est4; //4th estimate improvement for dxdt, dydt, and dzdt
float dydt_est4; //
float dzdt_est4; //
float simpsons; //variable used for Simpson's Rule
FILE *infile; //brings in parameter file, in.txt
FILE *outfile; //outputs results to file, out.txt
FILE *outfile2; //outputs results for y integral to out2.txt
infile = fopen("in.txt", "r"); // file input using relative path name of file
outfile = fopen("out.txt", "w"); // file output using relative path name of file
outfile2 = fopen("out2.txt", "w"); // file output using relative path name of file
fscanf(infile,"%f %f %f %f %f %f %f %f", &lambda, &mu, &x0, &y0, &z0, &tinit, &tfinal, &tdelta);
//scans paramemter values from input file and stores them to their repective variables
for (n = 0; n < 150; n++) //for loop to return values for x, y, and z at all 150 steps
{
dxdt_est = dxdt(lambda, x0, y0); //1st improvement using dxdt function from below
x_new = x0 + (tdelta * dxdt_est); //1st x value returned
dydt_est = dydt(lambda, x0, y0, mu);//1st improvement using dydt function from below
y_new = y0 + (tdelta * dydt_est); //1st y value returned
dzdt_est = dzdt(mu, y0); //1st improvement using dzdt function from below
z_new = z0 + (tdelta * dzdt_est); //1st z value returned
dxdt_est2 = ((dxdt_est + dxdt(lambda, x_new, y_new))/2); //2nd improvement using dxdt function from below
x_new2 = x_new + (tdelta * dxdt_est2); //2nd x value returned
dydt_est2 = ((dydt_est + dydt(lambda, x_new, y_new, mu))/2);//2nd improvement using dydt function from below
y_new2 = y_new + (tdelta * dydt_est2); //2nd y value returned
dzdt_est2 = ((dzdt_est + dzdt(mu, y_new))/2); //2nd improvement using dzdt function from below
z_new2 = z_new + (tdelta * dzdt_est2); //2nd z value returned
dxdt_est3 = ((dxdt_est2 + dxdt(lambda, x_new2, y_new2))/2); //3rd improvement using dxdt function from below
x_new3 = x_new2 + (tdelta * dxdt_est3); //3rd x value returned
dydt_est3 = ((dydt_est2 + dydt(lambda, x_new2, y_new2, mu))/2);//3rd improvement using dydt function from below
y_new3 = y_new2 + (tdelta * dydt_est3); //3rd y value returned
dzdt_est3 = ((dzdt_est2 + dzdt(mu, y_new2))/2); //3rd improvement using dzdt function from below
z_new3 = z_new2 + (tdelta * dzdt_est3); //3rd z value returned
dxdt_est4 = ((dxdt_est3 + dxdt(lambda, x_new3, y_new3))/2); //4th improvement using dxdt function from below
x_new4 = x_new3 + (tdelta * dxdt_est4); //4th x value returned
dydt_est4 = ((dydt_est3 + dydt(lambda, x_new3, y_new3, mu))/2);//4th improvement using dydt function from below
y_new4 = y_new3 + (tdelta * dydt_est4); //4th y value returned
dzdt_est4 = ((dzdt_est3 + dzdt(mu, y_new3))/2); //4th improvement using dzdt function from below
z_new4 = z_new3 + (tdelta * dzdt_est4); //4th z value returned
t_new = tinit + tdelta; //equation to step up time by 0.1 days every loop
simpsons = simp_y(t_new, y_new4); //calls function to give integral of y using Simpson's Rule
fprintf(outfile2, "%f\t %f\n", simpsons, t_new); //stores integral values of y to out2.txt file
float solution[3] = {x_new4, y_new4, z_new4}; //array definition to store x, y, and z values
printf("%f \t%f \t%f\n", solution[0], solution[1], solution[2]); //prints values for x_new4, y_new4, and z_new4
fprintf(outfile, "%f \t %f \t %f \t %f\n", solution[0], solution[1], solution[2], dydt_est4);
//outputs values of x_new4, y_new4, z_new4, and dydt_est4 into output file, out.txt
tinit = t_new; //equations to reassign values and continue loop
x0 = x_new;
y0 = y_new;
z0 = z_new;
x_new = x_new2;
y_new = y_new2;
z_new = z_new2;
x_new2 = x_new3;
y_new2 = y_new3;
z_new2 = z_new3;
x_new3 = x_new4;
y_new3 = y_new4;
z_new3 = z_new4;
};
}
float dxdt (float lambda, float y_new, float x_new) //dxdt base function
{
return -1 * lambda * y_new * x_new;
}
float dydt (float lambda ,float x_new, float y_new, float mu) //dydt base function
{
return (lambda * x_new * y_new) - (mu * y_new);
}
float dzdt(float mu,float y_new) //dzdt base function
{
return mu * y_new;
}
float simp_y (float t_new, float y_new4) //Function used to find integral of y using Simpson's Rule
{
return (((t_new-(t_new-2))/2) * (((y_new4)-2) + (4 * ((y_new4 + (y_new4 - 2))/2)) + y_new4));
}
"The instructor gave me this guidance to fix some things:
