#include <numeric.h> #include <stdio.h> #include <chplot.h> #define NVAR 4 #define POINTS 50 int main() { double t0=1, tf=10, y0[NVAR]; void func(double t, double y[], double dydt[], void *param); double t[POINTS], y[NVAR][POINTS]; int i, points; class CPlot plot; y0[0]=j0(t0); y0[1]=j1(t0); y0[2]=jn(2,t0); y0[3]=jn(3,t0); points = oderk(func, t0, tf, y0, NULL, t, y); printf("points = %d\n", points); printf("\n%8s %18s %15s\n","t","integral","bessj(2,t)"); for (i=0; i<points; i++) printf("%10.4f %14.6f %14.6f\n", t[i],y[2][i],jn(2,t[i])); printf("\nThe leftover space in t and y are padded with results at tf\n"); for (i=points; i<POINTS; i++) printf("%10.4f %14.6f %14.6f\n", t[i],y[2][i],jn(2,t[i])); plotxy(t, y, "solving ode using oderk()", "t", "yi", &plot); plot.legend("j0", 0); plot.legend("j1", 1); plot.legend("j2", 2); plot.legend("j3", 3); plot.plotting(); } void func(double t, double y[], double dydt[], void *param) { dydt[0] = -y[1]; dydt[1]=y[0]-(1.0/t)*y[1]; dydt[2]=y[1]-(2.0/t)*y[2]; dydt[3]=y[2]-(3.0/t)*y[3]; }
points = 43 t integral bessj(2,t) 1.0000 0.114903 0.114903 1.1995 0.159225 0.159225 1.4208 0.212473 0.212473 1.6107 0.259628 0.259628 1.8230 0.311687 0.311687 2.0595 0.365944 0.365944 2.3217 0.417776 0.417776 2.5117 0.447689 0.447689 2.7100 0.470487 0.470487 2.9266 0.484255 0.484255 3.1546 0.485094 0.485094 3.3763 0.471942 0.471942 3.5876 0.446661 0.446661 3.7842 0.412467 0.412467 4.0366 0.354910 0.354910 4.2242 0.303546 0.303546 4.4402 0.237267 0.237267 4.6658 0.162285 0.162285 4.8797 0.088340 0.088340 5.0780 0.019679 0.019679 5.3301 -0.064486 -0.064486 5.5121 -0.120926 -0.120926 5.7491 -0.186326 -0.186326 5.9414 -0.230973 -0.230973 6.1517 -0.269627 -0.269627 6.3458 -0.294896 -0.294896 6.5806 -0.311270 -0.311270 6.7618 -0.313089 -0.313089 6.9592 -0.304516 -0.304516 7.2117 -0.278417 -0.278417 7.3951 -0.249841 -0.249841 7.6042 -0.208797 -0.208797 7.8362 -0.154879 -0.154879 8.0631 -0.096271 -0.096271 8.2728 -0.039407 -0.039407 8.5323 0.030949 0.030949 8.7792 0.094212 0.094212 9.0451 0.154300 0.154300 9.2558 0.193600 0.193600 9.4769 0.225200 0.225200 9.6799 0.244472 0.244472 9.9231 0.254481 0.254481 10.0000 0.254630 0.254630 The leftover space in t and y are padded with results at tf 10.0000 0.254630 0.254630 10.0000 0.254630 0.254630 10.0000 0.254630 0.254630 10.0000 0.254630 0.254630 10.0000 0.254630 0.254630 10.0000 0.254630 0.254630 10.0000 0.254630 0.254630