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 |
Assignment 2 - Question1 April 12, 2007 12:09PM |
Registered: 17 years ago Posts: 24 Rating: 0 |
x y(x) Anal Error 0 3 3 0 0.2 3.2213563563 3.2184741262 -0.0028822301292 0.4 3.478586502 3.4353877658 -0.043198736218 0.6 3.6756088018 3.4936904666 -0.18191833519 0.8 3.5185266651 3.1484556902 -0.37007097492 1 2.5937700615 2.4287508015 -0.16501926002 1.2 2.0123757724 2.0073658525 -0.0050099199651
x y(x) Anal Error 0 3 3 0 0.2 3.2219411639 3.2184741262 -0.003467037713 0.4 3.4827684135 3.4353877658 -0.047380647785 0.6 3.6934705102 3.4936904666 -0.19978004362 0.8 3.5615923979 3.1484556902 -0.41313670766 1 2.5536509554 2.4287508015 -0.12490015394 1.2 2.9159785505 2.0073658525 -0.90861269806
Re: Assignment 2 - Question1 April 14, 2007 11:42AM |
Registered: 17 years ago Posts: 2,813 Rating: 0 |
euler, x0 = 0, y0 = 3: iteration 2: x = 0.200, y = 3.2093399780, error = 0.0091341482 iteration 4: x = 0.400, y = 3.4353981802, error = 0.0000104144 iteration 6: x = 0.600, y = 3.5589858030, error = 0.0652953364 iteration 8: x = 0.800, y = 3.3171487422, error = 0.1686930520 iteration 10: x = 1.000, y = 2.5211080221, error = 0.0923572206 iteration 12: x = 1.200, y = 1.9439908411, error = 0.0633750113 taylor2, x0 = 0, y0 = 3: iteration 2: x = 0.200, y = 3.2193012972, error = 0.0008271710 iteration 4: x = 0.400, y = 3.4422398293, error = 0.0068520635 iteration 6: x = 0.600, y = 3.5121725352, error = 0.0184820686 iteration 8: x = 0.800, y = 3.1660699391, error = 0.0176142489 iteration 10: x = 1.000, y = 2.4110724759, error = 0.0176783256 iteration 12: x = 1.200, y = 2.0052309742, error = 0.0021348783 rk2, x0 = 0, y0 = 3: iteration 2: x = 0.200, y = 3.2180246082, error = 0.0004495181 iteration 4: x = 0.400, y = 3.4343646586, error = 0.0010231072 iteration 6: x = 0.600, y = 3.4912198381, error = 0.0024706285 iteration 8: x = 0.800, y = 3.1391799347, error = 0.0092757556 iteration 10: x = 1.000, y = 2.4227839691, error = 0.0059668323 iteration 12: x = 1.200, y = 2.1634986912, error = 0.1561328387 rk4, x0 = 0, y0 = 3: iteration 2: x = 0.200, y = 3.2184738615, error = 0.0000002647 iteration 4: x = 0.400, y = 3.4353864230, error = 0.0000013428 iteration 6: x = 0.600, y = 3.4936848385, error = 0.0000056281 iteration 8: x = 0.800, y = 3.1484295110, error = 0.0000261793 iteration 10: x = 1.000, y = 2.4290026395, error = 0.0002518381 iteration 12: x = 1.200, y = 2.1110075716, error = 0.1036417191
Re: Assignment 2 - Question1 April 15, 2007 01:21AM |
Registered: 17 years ago Posts: 24 Rating: 0 |
Taylor2 x y(x) Anal Error 0 3 3 0 0.2 3.2203434856 3.2214018652 0.0010583796276 0.4 3.4789340058 3.4918072474 0.012873241537 0.6 3.7619613351 3.8220109013 0.060049566221 0.8 4.0401487438 4.2251244334 0.18497568958 1 4.2662939357 4.7170399951 0.45074605948 1.2 4.3796849357 5.3169723305 0.93728739484 Taylor4 x y(x) Anal Error 0 3 3 0 0.2 3.2195963188 3.2214018652 0.0018055463819 0.4 3.4726316334 3.4918072474 0.019175614009 0.6 3.7413238168 3.8220109013 0.08068708453 0.8 3.9923358575 4.2251244334 0.23278857587 1 4.177000614 4.7170399951 0.54003938114 1.2 4.2396836985 5.3169723305 1.077288632 Euler x y(x) Anal Error 0 3 3 0 0.2 3.2099884808 3.2214018652 0.011413384388 0.4 3.4635982578 3.4918072474 0.028208989587 0.6 3.7677687253 3.8220109013 0.054242176011 0.8 4.1276312448 4.2251244334 0.097493188594 1 4.5441040762 4.7170399951 0.17293591897 1.2 5.0104163598 5.3169723305 0.30655597073 Modified Euler x y(x) Anal Error 0 3 3 0 0.2 3.2206867238 3.2214018652 0.00071514141174 0.4 3.4888235676 3.4918072474 0.0029836797582 0.6 3.8120863868 3.8220109013 0.0099245144939 0.8 4.1955071376 4.2251244334 0.029617295735 1 4.6387599313 4.7170399951 0.07828006386 1.2 5.1320192559 5.3169723305 0.18495307459 RK2 x y(x) Anal Error 0 3 3 0 0.2 3.2209324296 3.2214018652 0.00046943560882 0.4 3.4897092955 3.4918072474 0.0020979519022 0.6 3.8137307444 3.8220109013 0.0082801568964 0.8 4.1979304957 4.2251244334 0.027193937719 1 4.6417910729 4.7170399951 0.075248922261 1.2 5.1351920505 5.3169723305 0.18178028005 RK4 x y(x) Anal Error 0 3 3 0 0.2 3.2213514096 3.2214018652 5.0455588525e-005 0.4 3.4908247448 3.4918072474 0.00098250256804 0.6 3.8159461128 3.8220109013 0.006064788528 0.8 4.2018015061 4.2251244334 0.023322927313 1 4.6480354929 4.7170399951 0.069004502217 1.2 5.1446631998 5.3169723305 0.17230913073 RKF x y(x) Anal Error 0 3 3 0 0.2 3.2213515972 3.2214018652 5.0267986667e-005 0.4 3.4908252008 3.4918072474 0.00098204658103 0.6 3.8159469326 3.8220109013 0.0060639687363 0.8 4.2018027884 4.2251244334 0.023321644964 1 4.6480373209 4.7170399951 0.069002674257 1.2 5.1446656186 5.3169723305 0.17230671194
Taylor2 x y(x) Anal Error 0 3 3 0 0.2 3.2193012972 3.2184741262 -0.0008271710359 0.4 3.4422398293 3.4353877658 -0.0068520635245 0.6 3.5121725352 3.4936904666 -0.018482068568 0.8 3.1660699391 3.1484556902 -0.017614248919 1 2.4110724759 2.4287508015 0.017678325625 1.2 2.0052309742 2.0073658525 0.0021348782826 Taylor4 x y(x) Anal Error 0 3 3 0 0.2 3.2185207946 3.2184741262 -4.6668434721e-005 0.4 3.4358441109 3.4353877658 -0.00045634513064 0.6 3.4971656277 3.4936904666 -0.0034751610501 0.8 3.1598388275 3.1484556902 -0.011383137278 1 2.4330536317 2.4287508015 -0.0043028301897 1.2 1.9857634462 2.0073658525 0.021602406235 Euler x y(x) Anal Error 0 3 3 0 0.2 3.209339978 3.2184741262 0.0091341482124 0.4 3.4353981802 3.4353877658 -1.0414440053e-005 0.6 3.558985803 3.4936904666 -0.065295336361 0.8 3.3171487422 3.1484556902 -0.16869305196 1 2.5211080221 2.4287508015 -0.0923572206 1.2 1.9439908411 2.0073658525 0.063375011348 Modified Euler x y(x) Anal Error 0 3 3 0 0.2 3.2171252975 3.2184741262 0.0013488287112 0.4 3.4307174537 3.4353877658 0.0046703120859 0.6 3.4833159116 3.4936904666 0.010374555026 0.8 3.1257658241 3.1484556902 0.022689866163 1 2.4235454081 2.4287508015 0.0052053933709 1.2 2.0216064616 2.0073658525 -0.014240609186 RK2 x y(x) Anal Error 0 3 3 0 0.2 3.215721898 3.2184741262 0.0027522282176 0.4 3.4235160934 3.4353877658 0.011871672391 0.6 3.4647292163 3.4936904666 0.028961250299 0.8 3.1063473583 3.1484556902 0.042108331911 1 2.4454673608 2.4287508015 -0.016716559362 1.2 15.084976093 2.0073658525 -13.07761024 RK4 x y(x) Anal Error 0 3 3 0 0.2 3.2184738615 3.2184741262 2.6474063269e-007 0.4 3.435386423 3.4353877658 1.3427904195e-006 0.6 3.4936848385 3.4936904666 5.6281019506e-006 0.8 3.148429511 3.1484556902 2.617927303e-005 1 2.4290026395 2.4287508015 -0.00025183805774 1.2 2.1110075716 2.0073658525 -0.1036417191 RKF x y(x) Anal Error 0 3 3 0 0.2 3.2184741414 3.2184741262 -1.5206708621e-008 0.4 3.4353878194 3.4353877658 -5.3637161181e-008 0.6 3.4936901309 3.4936904666 3.3571779348e-007 0.8 3.1484542082 3.1484556902 1.4820231605e-006 1 2.4287720091 2.4287508015 -2.1207582356e-005 1.2 1.9655328481 2.0073658525 0.041833004336
Re: Assignment 2 - Question1 April 15, 2007 01:36AM |
Registered: 17 years ago Posts: 2,813 Rating: 0 |
Re: Assignment 2 - Question1 April 16, 2007 07:50AM |
Registered: 17 years ago Posts: 2,813 Rating: 0 |
#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; }