Assignment 2 - Question1
Posted by murfinke
| Announcements | Last Post | |
|---|---|---|
 |
SoC Curricula | 09/30/2017 01:08PM |
|
Demarcation or scoping of examinations and assessment | 02/13/2017 07:59AM |
|
School of Computing Short Learning Programmes | 11/24/2014 08:37AM |
|
Unisa contact information | 07/28/2011 01:28PM |
|
April 14, 2007 11:42AM |
|
April 15, 2007 01:36AM |
murfinke Wrote:
-------------------------------------------------------
> This might seem like a silly question but should
> we be computing our values in degrees or radians.
> C++ automaticaly uses radians. All your answers
> are computed for radians with the Trig functions
> (tan,cos...). If I convert to degrees I don't get
> such erratic values
it'll definitely flatten the function out a lot, but that's solving a different problem. i'm pretty sure a nasty function was chosen (i haven't plotted it) over a nasty interval so as to allow the strengths of the respective methods to shine through. my taylor2 implementation is particularly efficient, and i wish t4 didn't have such massive expressions to evaluate.
at some stage i need to look at rkf and a-m too... first some sleep (serious lack thereof lately due to cos301)! i'll digest those results then
-------------------------------------------------------
> This might seem like a silly question but should
> we be computing our values in degrees or radians.
> C++ automaticaly uses radians. All your answers
> are computed for radians with the Trig functions
> (tan,cos...). If I convert to degrees I don't get
> such erratic values
it'll definitely flatten the function out a lot, but that's solving a different problem. i'm pretty sure a nasty function was chosen (i haven't plotted it) over a nasty interval so as to allow the strengths of the respective methods to shine through. my taylor2 implementation is particularly efficient, and i wish t4 didn't have such massive expressions to evaluate.
at some stage i need to look at rkf and a-m too... first some sleep (serious lack thereof lately due to cos301)! i'll digest those results then
|
April 16, 2007 07:50AM |
it's quite long, but here's some cleaned up code. i find stepsizes 0.01 and 0.001 to be much more illustrative than 0.1, so they are used here.
the order estimation results (min,avg,max):
taylor2: (1.98, 2.01, 2.06)
euler: (0.96, 1.00, 1.02)
euler2: (1.99, 2.02, 2.12)
rk2: (1.35, 1.91, 2.10)
rk4: (3.86, 4.02, 4.12)
rkf: (1.57, 3.53, 5.12)
code:
the order estimation results (min,avg,max):
taylor2: (1.98, 2.01, 2.06)
euler: (0.96, 1.00, 1.02)
euler2: (1.99, 2.02, 2.12)
rk2: (1.35, 1.91, 2.10)
rk4: (3.86, 4.02, 4.12)
rkf: (1.57, 3.53, 5.12)
code:
#include <math.h>
#include <stdio.h>
typedef double (*deriv_func)(double x, double y);
inline double deriv_1(double x, double y)
{
return (1.0 - 6.0 * x * tan(x * x)) * (y - 2.0);
}
inline double deriv_2(double x, double y)
{
const double
x2 = x * x,
tan_term = tan(x2),
sec_term = x / cos(x2),
big_term = 1.0 - 6.0 * x * tan_term;
return (y - 2.0) * (big_term * big_term - 6.0 * (2.0 * sec_term * sec_term + tan_term));
}
double deriv_q2(double x, double y)
{
return (3.0 * y - 54.0) / x - 8.0 * x + 42.0;
}
inline double F_3(const double x)
{
// cubing via 3 muls is much faster than pow, but less accurate
return 2.0 + exp(x) * pow(cos(x * x),3.0);
}
double taylor_2(const double x0, const double y0, const double h,
const int num_iters, const int print_iters,
double * y_out, double * error_out)
{
const double h2 = h * h * 0.5;
double y = y0;
for (int i = 0, j = 0; i <= num_iters; ++i)
{
const double
x_n = x0 + h * double(i), // to prevent cumulative error in x_n
x2 = x_n * x_n,
y_2_term = y - 2.0,
tan_term = tan(x2),
sec_term = x_n / cos(x2),
big_term = 1.0 - 6.0 * x_n * tan_term,
d_1 = y_2_term * big_term,
d_2 = y_2_term * (big_term * big_term - 6.0 * (2.0 * sec_term * sec_term + tan_term));
y += h * d_1 + h2 * d_2;
//y += h * deriv_1(x_n,y) + h2 * deriv_2(x_n,y); // computation is folded above
if (!((i + 1) % print_iters))
{
y_out[j] = y; error_out[j] = y - F_3(x_n + h);
printf("x = %.1f, y = %.10f, error = %.10f\n",x_n + h,y,error_out[j]);
++j;
}
}
return y;
}
double euler_1( const double x0, const double y0, const double h,
const int num_iters, const int print_iters,
double * y_out, double * error_out)
{
double y = y0;
for (int i = 0, j = 0; i <= num_iters; ++i)
{
const double x_n = x0 + h * double(i);
y += h * deriv_1(x_n,y);
if (!((i + 1) % print_iters))
{
y_out[j] = y; error_out[j] = y - F_3(x_n + h);
printf("x = %.1f, y = %.10f, error = %.10f\n",x_n + h,y,error_out[j]);
++j;
}
}
return y;
}
double euler_2( const double x0, const double y0, const double h,
const int num_iters, const int print_iters,
double * y_out, double * error_out)
{
double y = y0;
for (int i = 0, j = 0; i <= num_iters; ++i)
{
const double
x_n = x0 + h * double(i),
x_n1 = x0 + h * double(i + 1),
y_dash = deriv_1(x_n,y),
y_dash1 = deriv_1(x_n1,y + y_dash * h),
y_dash_final = (y_dash + y_dash1) * 0.5;
y += h * y_dash_final;
if (!((i + 1) % print_iters))
{
y_out[j] = y; error_out[j] = y - F_3(x_n + h);
printf("x = %.1f, y = %.10f, error = %.10f\n",x_n + h,y,error_out[j]);
++j;
}
}
return y;
}
double rk2(const double x0, const double y0, const double h,
const int num_iters, const int print_iters,
double * y_out, double * error_out)
{
// use special(ised?) weights given in the tutorial letter
const double weights[2] = { h / 3.0, h * 2.0 / 3.0 };
const double alphabeta = h * 3.0 / 4.0;
double y = y0;
for (int i = 0, j = 0; i <= num_iters; ++i)
{
const double
x_n = x0 + h * double(i),
k_1 = deriv_1(x_n,y),
k_2 = deriv_1(x_n + alphabeta,y + alphabeta * k_1);
y += k_1 * weights[0] + k_2 * weights[1];
if (!((i + 1) % print_iters))
{
y_out[j] = y; error_out[j] = y - F_3(x_n + h);
printf("x = %.1f, y = %.10f, error = %.10f\n",x_n + h,y,error_out[j]);
++j;
}
}
return y;
}
double rk4(const double x0, const double y0, const double h,
const int num_iters, const int print_iters,
double * y_out, double * error_out, deriv_func deriv_f)
{
const double consts[2] = { h / 2.0, h / 6.0 };
double y = y0;
for (int i = 0, j = 0; i <= num_iters; ++i)
{
const double
x_n = x0 + h * double(i),
x_2 = x_n + consts[0],
k_1 = deriv_f(x_n,y),
k_2 = deriv_f(x_2,y + k_1 * consts[0]),
k_3 = deriv_f(x_2,y + k_2 * consts[0]),
k_4 = deriv_f(x_n + h,y + k_3 * h);
y += (k_1 + k_2 * 2.0 + k_3 * 2.0 + k_4) * consts[1];
if (!((i + 1) % print_iters))
{
y_out[j] = y; error_out[j] = y - F_3(x_n + h);
printf("x = %.1f, y = %.10f, error = %.10f\n",x_n + h,y,error_out[j]);
++j;
}
}
return y;
}
double rkf(const double x0, const double y0, const double h,
const int num_iters, const int print_iters,
double * y_out, double * error_out)
{
// this mass of constants should be put in a more global scope if rkf will be called often,
// so they only need to be recomputed when the stepsize h changes
const double weights[] =
{
h / 4.0,
h * 3.0 / 8.0, h * 3.0 / 32.0, h * 9.0 / 32.0,
h * 12.0 / 13.0, h * 1932.0 / 2197.0, h * -7200.0 / 2197.0, h * 7296.0 / 2197.0,
h, h * 439.0 / 216.0, h * -8.0, h * 3680.0 / 513.0, h * -845.0 / 4104.0,
h / 2.0, h * -8.0 / 27.0, h * 2.0, h * -3544.0 / 2565.0, h * 1859.0 / 4104.0, h * -11.0 / 40.0,
h * 25.0 / 216.0, h * 1408.0 / 2565.0, h * 2197.0 / 4104.0, h / -5.0,
h * 16.0 / 135.0, h * 6656.0 / 12825.0, h * 28561.0 / 56430.0, h * -9.0 / 50.0, h * 2.0 / 55.0,
h / 360.0, h * -128.0 / 4275.0, h * -2197.0 / 75420.0, h / 50.0, h * 2.0 / 55.0
};
double y = y0;
for (int i = 0, j = 0; i <= num_iters; ++i)
{
const double
x_n = x0 + h * double(i),
k_1 = deriv_1(x_n,y),
k_2 = deriv_1(x_n + weights[ 0],y + weights[ 0] * k_1),
k_3 = deriv_1(x_n + weights[ 1],y + weights[ 2] * k_1 + weights[ 3] * k_2),
k_4 = deriv_1(x_n + weights[ 4],y + weights[ 5] * k_1 + weights[ 6] * k_2 + weights[ 7] * k_3),
k_5 = deriv_1(x_n + weights[ 8],y + weights[ 9] * k_1 + weights[10] * k_2 + weights[11] * k_3 +
weights[12] * k_4),
k_6 = deriv_1(x_n + weights[13],y + weights[14] * k_1 + weights[15] * k_2 + weights[16] * k_3 +
weights[17] * k_4 + weights[18] * k_5);
// even if it's not used, might as well show its computation
const double y_hat = y + weights[19] * k_1 + weights[20] * k_3 + weights[21] * k_4 + weights[22] * k_5;
y += weights[23] * k_1 + weights[24] * k_3 + weights[25] * k_4 + weights[26] * k_5 + weights[27] * k_6;
if (!((i + 1) % print_iters))
{
y_out[j] = y; error_out[j] = y - F_3(x_n + h);
printf("x = %.1f, y = %.10f, error = %.10f\n",x_n + h,y,error_out[j]);
++j;
}
}
return y;
}
// utility function for gathering stats
void minavgmax(const int howmany, double *data_in, double *data_out)
{
double minval = -log(0.0), maxval = log(0.0), avgval = 0.0;
for (int i = 0; i < howmany; ++i)
{
avgval += data_in;
minval = (minval < data_in) ? minval : data_in;
maxval = (maxval > data_in) ? maxval : data_in;
}
data_out[0] = minval;
data_out[1] = avgval / double(howmany);
data_out[2] = maxval;
}
// should make all this dynamic, likewise the derivative function
#define METHODS 6
#define RESULTS 6
#define STEPSIZES 2
int main(int argc, char ** argv)
{
int i, j;
double step_sizes[STEPSIZES] = { 0.01, 0.001 };
double *y_out = new double[METHODS*RESULTS*STEPSIZES];
double *e_out = new double[METHODS*RESULTS*STEPSIZES];
for (i = 0; i < STEPSIZES; ++i)
{
printf("using stepsize = %.3f\n\n",step_sizes);
const double
x0 = 0.0,
y0 = 3.0,
x_interval = 1.2, // solution interval
p_interval = 0.2; // printing interval
const int
num_iters = ceil(x_interval / step_sizes),
prt_iters = ceil(p_interval / step_sizes);
double
*y_batch = &y_out[i * METHODS * RESULTS],
*e_batch = &e_out[i * METHODS * RESULTS],
step = step_sizes;
printf("taylor expansion, 2nd order:\n"
;
taylor_2(x0,y0,step,num_iters,prt_iters,&y_batch[0*METHODS],&e_batch[0*METHODS]);
printf("\n";
printf("euler, 1st order:\n";
euler_1(x0,y0,step,num_iters,prt_iters,&y_batch[1*METHODS],&e_batch[1*METHODS]);
printf("\n";
printf("modified euler, 2nd order:\n";
euler_2(x0,y0,step,num_iters,prt_iters,&y_batch[2*METHODS],&e_batch[2*METHODS]);
printf("\n";
printf("runge-kutta, 2nd order:\n";
rk2(x0,y0,step,num_iters,prt_iters,&y_batch[3*METHODS],&e_batch[3*METHODS]);
printf("\n";
printf("runge-kutta, 4th order:\n";
rk4(x0,y0,step,num_iters,prt_iters,&y_batch[4*METHODS],&e_batch[4*METHODS],deriv_1);
printf("\n";
printf("runge-kutta-fehlberg, 5th order:\n";
rkf(x0,y0,step,num_iters,prt_iters,&y_batch[5*METHODS],&e_batch[5*METHODS]);
printf("\n";
}
const double log_scale = log(step_sizes[0] / step_sizes[1]);
double orders[METHODS * RESULTS], stats[3];
for (j = 0; j < METHODS; ++j)
for (i = 0; i < RESULTS; ++i)
orders[j * RESULTS + i] = log(abs(e_out[j * RESULTS + i]) /
abs(e_out[j * RESULTS + i + METHODS*RESULTS])) / log_scale;
printf("estimates for the error order (min,avg,max):\n\n";
minavgmax(RESULTS,&orders[0*METHODS],(double *)&stats[0]);
printf("taylor2: (%.2f, %.2f, %.2f)\n",stats[0],stats[1],stats[2]);
minavgmax(RESULTS,&orders[1*METHODS],(double *)&stats);
printf("euler: (%.2f, %.2f, %.2f)\n",stats[0],stats[1],stats[2]);
minavgmax(RESULTS,&orders[2*METHODS],(double *)&stats);
printf("euler2: (%.2f, %.2f, %.2f)\n",stats[0],stats[1],stats[2]);
minavgmax(RESULTS,&orders[3*METHODS],(double *)&stats);
printf("rk2: (%.2f, %.2f, %.2f)\n",stats[0],stats[1],stats[2]);
minavgmax(RESULTS,&orders[4*METHODS],(double *)&stats);
printf("rk4: (%.2f, %.2f, %.2f)\n",stats[0],stats[1],stats[2]);
minavgmax(RESULTS,&orders[5*METHODS],(double *)&stats);
printf("rkf: (%.2f, %.2f, %.2f)\n",stats[0],stats[1],stats[2]);
// now gather min/avg/max stats for all stepsize-pairs...
delete [] y_out;
delete [] e_out;
return 0;
}
Sorry, only registered users may post in this forum.