「プログラムソース」
７−１、-------------------------------------------------------------
#include <stdio.h>
#include <math.h>

#define MAX 100

double f1(double t,double x_1,double x_2,double x_3,double x_4);
double f2(double t,double x_1,double x_2,double x_3,double x_4);
double f3(double t,double x_1,double x_2,double x_3,double x_4);
double f4(double t,double x_1,double x_2,double x_3,double x_4);

int main(void){

		double t,dt;
		double x_1,x_2,x_3,x_4;
		double k0[4],k1[4],k2[4],k3[4];

		FILE *rkdata;
		rkdata=fopen("rkdata.txt","w");
		
		dt=0.01;
		x_1=1;
		x_2=0;
		x_3=0;
		x_4=1;


		printf("rkdata.txtを実行ファイルと同じフォルダに出力しました。\n");

		for(t=0.0;t<MAX;t+=dt) {
			k0[0]=dt*f1(t,x_1,x_2,x_3,x_4);
			k0[1]=dt*f2(t,x_1,x_2,x_3,x_4);
			k0[2]=dt*f3(t,x_1,x_2,x_3,x_4);
			k0[3]=dt*f4(t,x_1,x_2,x_3,x_4);

			k1[0]=dt*f1(t+dt/2.0,x_1+k0[0]/2.0,x_2+k0[1]/2.0,x_3+k0[2]/2.0,x_4+k0[3]/2.0);
			k1[1]=dt*f2(t+dt/2.0,x_1+k0[0]/2.0,x_2+k0[1]/2.0,x_3+k0[2]/2.0,x_4+k0[3]/2.0);
			k1[2]=dt*f3(t+dt/2.0,x_1+k0[0]/2.0,x_2+k0[1]/2.0,x_3+k0[2]/2.0,x_4+k0[3]/2.0);
			k1[3]=dt*f4(t+dt/2.0,x_1+k0[0]/2.0,x_2+k0[1]/2.0,x_3+k0[2]/2.0,x_4+k0[3]/2.0);
				
			k2[0]=dt*f1(t+dt/2.0,x_1+k1[0]/2.0,x_2+k1[1]/2.0,x_3+k1[2]/2.0,x_4+k1[3]/2.0);
			k2[1]=dt*f2(t+dt/2.0,x_1+k1[0]/2.0,x_2+k1[1]/2.0,x_3+k1[2]/2.0,x_4+k1[3]/2.0);
			k2[2]=dt*f3(t+dt/2.0,x_1+k1[0]/2.0,x_2+k1[1]/2.0,x_3+k1[2]/2.0,x_4+k1[3]/2.0);
			k2[3]=dt*f4(t+dt/2.0,x_1+k1[0]/2.0,x_2+k1[1]/2.0,x_3+k1[2]/2.0,x_4+k1[3]/2.0);
				
			k3[0]=dt*f1(t+dt,x_1+k2[0],x_2+k2[1],x_3+k2[2],x_4+k2[3]);
			k3[1]=dt*f2(t+dt,x_1+k2[0],x_2+k2[1],x_3+k2[2],x_4+k2[3]);
			k3[2]=dt*f3(t+dt,x_1+k2[0],x_2+k2[1],x_3+k2[2],x_4+k2[3]);
			k3[3]=dt*f4(t+dt,x_1+k2[0],x_2+k2[1],x_3+k2[2],x_4+k2[3]);

			x_1=x_1+(k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
			x_2=x_2+(k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
			x_3=x_3+(k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;
			x_4=x_4+(k0[3]+2.0*k1[3]+2.0*k2[3]+k3[3])/6.0;
	            
			fprintf(rkdata,"x={%f,%f,%f,%f}\n",x_1,x_2,x_3,x_4);

		}

		fclose(rkdata);

		return 0;
}

double f1(double t,double x_1,double x_2,double x_3,double x_4)
{
		double y;

		y=x_3;

		return y;
}

double f2(double t,double x_1,double x_2,double x_3,double x_4)
{
		double y;

		y=x_4;

		return y;
}

double f3(double t,double x_1,double x_2,double x_3,double x_4)
{
		double y,r;

		r=sqrt(x_1*x_1+x_2*x_2);
		y=-x_1/r*r*r;

		return y;
}

double f4(double t,double x_1,double x_2,double x_3,double x_4)
{
		double y,r;

		r=sqrt(x_1*x_1+x_2*x_2);
		y=-x_1/r*r*r;

		return y;
}



--------------------------------------------------------------------------------

「実行結果」
７−１--------------------------------------------------------------------------

実行結果は、膨大な量になったので、実行する事でテキストファイルに別途出力形式にしておきました。
テキストファイルに出力されてデータを一部コピペで載せておきます。


C:\Documents and Settings\i865gpen4\My Documents\プログラムソース>runngekutta
rkdata.txtを実行ファイルと同じフォルダに出力しました。
 -- Press any key to exit (Input "c" to continue) --

以下rkdata.txtの内容の一部をコピペ↓

x={0.999950,0.009950,-0.010000,0.990000}
x={0.999800,0.019800,-0.019999,0.980001}
x={0.999550,0.029550,-0.029995,0.970005}
x={0.999200,0.039200,-0.039989,0.960011}
x={0.998750,0.048750,-0.049978,0.950022}
x={0.998201,0.058201,-0.059962,0.940038}
x={0.997551,0.067551,-0.069940,0.930060}
x={0.996802,0.076802,-0.079910,0.920090}
x={0.995953,0.085953,-0.089871,0.910129}
x={0.995004,0.095004,-0.099822,0.900178}
x={0.993956,0.103956,-0.109761,0.890239}
x={0.992809,0.112809,-0.119688,0.880312}
x={0.991563,0.121563,-0.129601,0.870399}
x={0.990217,0.130217,-0.139499,0.860501}
x={0.988773,0.138773,-0.149380,0.850620}
x={0.987230,0.147230,-0.159243,0.840757}
x={0.985588,0.155588,-0.169087,0.830913}
x={0.983848,0.163848,-0.178911,0.821089}
x={0.982010,0.172010,-0.188713,0.811287}
x={0.980074,0.180074,-0.198491,0.801509}
x={0.978040,0.188040,-0.208245,0.791755}
x={0.975909,0.195909,-0.217972,0.782028}
x={0.973681,0.203681,-0.227672,0.772328}
x={0.971356,0.211356,-0.237343,0.762657}
x={0.968934,0.218934,-0.246984,0.753016}
x={0.966416,0.226416,-0.256593,0.743407}
x={0.963802,0.233802,-0.266168,0.733832}
x={0.961093,0.241093,-0.275709,0.724291}
x={0.958288,0.248288,-0.285214,0.714786}
x={0.955389,0.255389,-0.294682,0.705318}
x={0.952395,0.262395,-0.304110,0.695890}
x={0.949307,0.269307,-0.313498,0.686502}
x={0.946125,0.276125,-0.322845,0.677155}
x={0.942850,0.282850,-0.332148,0.667852}
x={0.939482,0.289482,-0.341406,0.658594}
x={0.936022,0.296022,-0.350619,0.649381}
x={0.932470,0.302470,-0.359784,0.640216}
x={0.928826,0.308826,-0.368900,0.631100}
x={0.925092,0.315092,-0.377967,0.622033}
x={0.921267,0.321267,-0.386981,0.613019}
x={0.917353,0.327353,-0.395943,0.604057}
x={0.913349,0.333349,-0.404851,0.595149}
x={0.909256,0.339256,-0.413704,0.586296}
x={0.905075,0.345075,-0.422499,0.577501}
x={0.900806,0.350806,-0.431237,0.568763}
x={0.896450,0.356450,-0.439915,0.560085}
x={0.892008,0.362008,-0.448533,0.551467}
x={0.887480,0.367480,-0.457089,0.542911}
x={0.882866,0.372866,-0.465582,0.534418}
x={0.878168,0.378168,-0.474011,0.525989}
x={0.873386,0.383386,-0.482374,0.517626}
x={0.868521,0.388521,-0.490672,0.509328}
x={0.863573,0.393573,-0.498901,0.501099}
x={0.858543,0.398543,-0.507062,0.492938}
x={0.853432,0.403432,-0.515154,0.484846}
x={0.848240,0.408240,-0.523174,0.476826}
x={0.842969,0.412969,-0.531124,0.468876}
x={0.837618,0.417618,-0.539000,0.461000}
x={0.832189,0.422189,-0.546803,0.453197}
x={0.826682,0.426682,-0.554531,0.445469}
x={0.821099,0.431099,-0.562184,0.437816}
x={0.815439,0.435439,-0.569760,0.430240}
x={0.809704,0.439704,-0.577259,0.422741}
x={0.803894,0.443894,-0.584681,0.415319}
x={0.798010,0.448010,-0.592024,0.407976}
x={0.792054,0.452054,-0.599287,0.400713}
x={0.786025,0.456025,-0.606470,0.393530}
x={0.779925,0.459925,-0.613572,0.386428}
x={0.773754,0.463754,-0.620593,0.379407}
x={0.767513,0.467513,-0.627532,0.372468}
x={0.761203,0.471203,-0.634388,0.365612}
x={0.754826,0.474826,-0.641161,0.358839}
x={0.748380,0.478380,-0.647850,0.352150}
x={0.741869,0.481869,-0.654455,0.345545}
x={0.735292,0.485292,-0.660976,0.339024}
x={0.728650,0.488650,-0.667411,0.332589}
x={0.721944,0.491944,-0.673761,0.326239}
x={0.715175,0.495175,-0.680025,0.319975}
x={0.708343,0.498343,-0.686203,0.313797}
x={0.701451,0.501451,-0.692295,0.307705}
x={0.694498,0.504498,-0.698300,0.301700}
x={0.687485,0.507485,-0.704218,0.295782}
x={0.680414,0.510414,-0.710049,0.289951}
x={0.673285,0.513285,-0.715793,0.284207}
x={0.666098,0.516098,-0.721449,0.278551}
x={0.658856,0.518856,-0.727018,0.272982}
x={0.651558,0.521558,-0.732500,0.267500}
x={0.644206,0.524206,-0.737894,0.262106}
x={0.636801,0.526801,-0.743201,0.256799}
x={0.629342,0.529342,-0.748420,0.251580}
x={0.621832,0.531832,-0.753552,0.246448}
x={0.614272,0.534272,-0.758596,0.241404}
x={0.606661,0.536661,-0.763553,0.236447}
x={0.599001,0.539001,-0.768424,0.231576}
x={0.591293,0.541293,-0.773207,0.226793}
x={0.583537,0.543537,-0.777904,0.222096}
x={0.575735,0.545735,-0.782514,0.217486}
x={0.567887,0.547887,-0.787038,0.212962}
x={0.559994,0.549994,-0.791476,0.208524}
x={0.552058,0.552058,-0.795829,0.204171}
x={0.544078,0.554078,-0.800097,0.199903}
x={0.536056,0.556056,-0.804279,0.195721}
x={0.527993,0.557993,-0.808378,0.191622}
x={0.519889,0.559889,-0.812392,0.187608}
x={0.511745,0.561745,-0.816322,0.183678}
x={0.503563,0.563563,-0.820169,0.179831}
x={0.495342,0.565342,-0.823934,0.176066}
x={0.487084,0.567084,-0.827616,0.172384}
x={0.478790,0.568790,-0.831216,0.168784}
x={0.470460,0.570460,-0.834735,0.165265}
x={0.462096,0.572096,-0.838174,0.161826}
x={0.453697,0.573697,-0.841532,0.158468}
x={0.445265,0.575265,-0.844811,0.155189}
x={0.436801,0.576801,-0.848011,0.151989}
x={0.428305,0.578305,-0.851132,0.148868}
x={0.419779,0.579779,-0.854175,0.145825}
x={0.411222,0.581222,-0.857142,0.142858}
x={0.402636,0.582636,-0.860031,0.139969}
x={0.394022,0.584022,-0.862845,0.137155}
x={0.385379,0.585379,-0.865583,0.134417}
x={0.376710,0.586710,-0.868247,0.131753}
x={0.368015,0.588015,-0.870837,0.129163}
x={0.359294,0.589294,-0.873353,0.126647}
x={0.350548,0.590548,-0.875796,0.124204}
x={0.341778,0.591778,-0.878168,0.121832}
x={0.332985,0.592985,-0.880468,0.119532}
x={0.324169,0.594169,-0.882697,0.117303}
x={0.315331,0.595331,-0.884856,0.115144}
x={0.306472,0.596472,-0.886946,0.113054}
x={0.297592,0.597592,-0.888967,0.111033}
x={0.288693,0.598693,-0.890920,0.109080}
x={0.279774,0.599774,-0.892805,0.107195}
x={0.270837,0.600837,-0.894623,0.105377}
x={0.261882,0.601882,-0.896375,0.103625}
x={0.252910,0.602910,-0.898061,0.101939}
x={0.243921,0.603921,-0.899682,0.100318}
x={0.234916,0.604916,-0.901239,0.098761}
x={0.225896,0.605896,-0.902731,0.097269}
x={0.216862,0.606862,-0.904161,0.095839}
x={0.207813,0.607813,-0.905527,0.094473}
x={0.198752,0.608752,-0.906831,0.093169}
x={0.189677,0.609677,-0.908072,0.091928}
x={0.180590,0.610590,-0.909253,0.090747}
x={0.171492,0.611492,-0.910372,0.089628}
x={0.162383,0.612383,-0.911431,0.088569}
x={0.153264,0.613264,-0.912430,0.087570}
x={0.144135,0.614135,-0.913369,0.086631}
x={0.134997,0.614997,-0.914249,0.085751}
x={0.125850,0.615850,-0.915069,0.084931}
x={0.116695,0.616695,-0.915831,0.084169}
x={0.107534,0.617534,-0.916534,0.083466}
x={0.098365,0.618365,-0.917179,0.082821}
x={0.089190,0.619190,-0.917766,0.082234}
x={0.080010,0.620010,-0.918295,0.081705}
x={0.070824,0.620824,-0.918766,0.081234}
x={0.061635,0.621635,-0.919180,0.080820}
x={0.052441,0.622441,-0.919536,0.080464}
x={0.043244,0.623244,-0.919835,0.080165}
x={0.034045,0.624045,-0.920077,0.079923}
x={0.024843,0.624843,-0.920261,0.079739}
x={0.015639,0.625639,-0.920387,0.079613}
x={0.006435,0.626435,-0.920456,0.079544}
x={-0.002769,0.627231,-0.920468,0.079532}
x={-0.011974,0.628026,-0.920422,0.079578}
x={-0.021178,0.628822,-0.920317,0.079683}
x={-0.030380,0.629620,-0.920155,0.079845}
x={-0.039581,0.630419,-0.919934,0.080066}
x={-0.048779,0.631221,-0.919655,0.080345}
x={-0.057973,0.632027,-0.919317,0.080683}
x={-0.067165,0.632835,-0.918919,0.081081}
x={-0.076352,0.633648,-0.918462,0.081538}
x={-0.085534,0.634466,-0.917944,0.082056}
x={-0.094710,0.635290,-0.917366,0.082634}
x={-0.103881,0.636119,-0.916727,0.083273}
x={-0.113045,0.636955,-0.916027,0.083973}
x={-0.122201,0.637799,-0.915265,0.084735}
x={-0.131350,0.638650,-0.914440,0.085560}
x={-0.140490,0.639510,-0.913552,0.086448}
x={-0.149621,0.640379,-0.912600,0.087400}
x={-0.158742,0.641258,-0.911583,0.088417}
x={-0.167852,0.642148,-0.910502,0.089498}
x={-0.176951,0.643049,-0.909355,0.090645}
x={-0.186039,0.643961,-0.908142,0.091858}
x={-0.195114,0.644886,-0.906861,0.093139}
x={-0.204176,0.645824,-0.905512,0.094488}
x={-0.213224,0.646776,-0.904095,0.095905}
x={-0.222258,0.647742,-0.902608,0.097392}
x={-0.231276,0.648724,-0.901050,0.098950}
x={-0.240278,0.649722,-0.899422,0.100578}
x={-0.249264,0.650736,-0.897721,0.102279}
x={-0.258233,0.651767,-0.895947,0.104053}
x={-0.267183,0.652817,-0.894100,0.105900}
x={-0.276114,0.653886,-0.892178,0.107822}
x={-0.285026,0.654974,-0.890180,0.109820}
x={-0.293918,0.656082,-0.888106,0.111894}
x={-0.302788,0.657212,-0.885954,0.114046}
x={-0.311636,0.658364,-0.883723,0.116277}
x={-0.320462,0.659538,-0.881414,0.118586}
x={-0.329264,0.660736,-0.879023,0.120977}
x={-0.338042,0.661958,-0.876552,0.123448}
x={-0.346795,0.663205,-0.873998,0.126002}
x={-0.355522,0.664478,-0.871361,0.128639}
x={-0.364222,0.665778,-0.868639,0.131361}
x={-0.372895,0.667105,-0.865832,0.134168}
x={-0.381539,0.668461,-0.862939,0.137061}
x={-0.390153,0.669847,-0.859959,0.140041}
x={-0.398737,0.671263,-0.856890,0.143110}
x={-0.407291,0.672709,-0.853732,0.146268}
x={-0.415812,0.674188,-0.850484,0.149516}
x={-0.424300,0.675700,-0.847144,0.152856}
x={-0.432754,0.677246,-0.843713,0.156287}
x={-0.441174,0.678826,-0.840188,0.159812}
x={-0.449558,0.680442,-0.836569,0.163431}
x={-0.457905,0.682095,-0.832855,0.167145}
x={-0.466215,0.683785,-0.829045,0.170955}
x={-0.474485,0.685515,-0.825138,0.174862}
x={-0.482717,0.687283,-0.821133,0.178867}
x={-0.490908,0.689092,-0.817029,0.182971}
x={-0.499057,0.690943,-0.812826,0.187174}
x={-0.507164,0.692836,-0.808522,0.191478}
x={-0.515227,0.694773,-0.804116,0.195884}
x={-0.523246,0.696754,-0.799608,0.200392}
x={-0.531219,0.698781,-0.794997,0.205003}
x={-0.539146,0.700854,-0.790282,0.209718}
x={-0.547024,0.702976,-0.785462,0.214538}
x={-0.554854,0.705146,-0.780537,0.219463}
x={-0.562635,0.707365,-0.775505,0.224495}
x={-0.570364,0.709636,-0.770366,0.229634}
x={-0.578042,0.711958,-0.765119,0.234881}
x={-0.585666,0.714334,-0.759764,0.240236}
x={-0.593237,0.716763,-0.754299,0.245701}
x={-0.600752,0.719248,-0.748724,0.251276}
x={-0.608211,0.721789,-0.743039,0.256961}
x={-0.615612,0.724388,-0.737242,0.262758}
x={-0.622955,0.727045,-0.731334,0.268666}
x={-0.630239,0.729761,-0.725314,0.274686}
x={-0.637461,0.732539,-0.719180,0.280820}
x={-0.644622,0.735378,-0.712933,0.287067}
x={-0.651719,0.738281,-0.706572,0.293428}
x={-0.658753,0.741247,-0.700097,0.299903}
x={-0.665721,0.744279,-0.693507,0.306493}
x={-0.672623,0.747377,-0.686801,0.313199}
x={-0.679457,0.750543,-0.679980,0.320020}
x={-0.686222,0.753778,-0.673044,0.326956}
x={-0.692917,0.757083,-0.665990,0.334010}
x={-0.699541,0.760459,-0.658821,0.341179}
x={-0.706093,0.763907,-0.651534,0.348466}
x={-0.712572,0.767428,-0.644130,0.355870}
x={-0.718975,0.771025,-0.636609,0.363391}
x={-0.725303,0.774697,-0.628971,0.371029}
x={-0.731554,0.778446,-0.621215,0.378785}
x={-0.737727,0.782273,-0.613341,0.386659}
x={-0.743821,0.786179,-0.605350,0.394650}
x={-0.749834,0.790166,-0.597241,0.402759}
x={-0.755765,0.794235,-0.589014,0.410986}
x={-0.761614,0.798386,-0.580669,0.419331}
x={-0.767378,0.802622,-0.572206,0.427794}
x={-0.773057,0.806943,-0.563626,0.436374}
x={-0.778650,0.811350,-0.554929,0.445071}
x={-0.784156,0.815844,-0.546114,0.453886}
x={-0.789572,0.820428,-0.537182,0.462818}
x={-0.794899,0.825101,-0.528133,0.471867}
x={-0.800134,0.829866,-0.518968,0.481032}
x={-0.805278,0.834722,-0.509686,0.490314}
x={-0.810328,0.839672,-0.500288,0.499712}
x={-0.815283,0.844717,-0.490774,0.509226}
x={-0.820143,0.849857,-0.481145,0.518855}
x={-0.824906,0.855094,-0.471402,0.528598}
x={-0.829571,0.860429,-0.461544,0.538456}
x={-0.834136,0.865864,-0.451572,0.548428}
x={-0.838602,0.871398,-0.441486,0.558514}
x={-0.842966,0.877034,-0.431288,0.568712}
x={-0.847227,0.882773,-0.420978,0.579022}
x={-0.851385,0.888615,-0.410556,0.589444}
x={-0.855438,0.894562,-0.400023,0.599977}
x={-0.859385,0.900615,-0.389380,0.610620}
x={-0.863225,0.906775,-0.378627,0.621373}
x={-0.866957,0.913043,-0.367766,0.632234}



--------------------------------------------------------------------------------

「考察」
オイラー法とは、微分方程式dx/dt=f(x,t)を数値的に解くために、
連続な微分を有限幅の差分にすることによって求めるやり方で、

dx/dt ~= (x(t+�冲)-x(t))/�冲

のように近似して問題を解く。まず、微分を時間ｔにおける
ものとして解釈すると以下のようになる。

(x(t+�冲)-x(t))/�冲 = f(x(t),t)

より、時間ｔのｘが与えられたときに時間(t+�冲)のｘは、

x(t+�冲)=x(t)+f(x(t),t)�冲

と計算することができる。すなわち、時間をｎ分割すると、

x[n+1]=x[n]+f(x[n],t)�冲

により、時間積分を行うことができる。

このようなオイラー法から始まり、ホイン法といったように、
パラメーターを増やして誤差項の次数を上げていく方法で、
最良の方法と言われるのが4次のルンゲクッタ法だとされている。

パラメーターを増やして、5, 6, 7, と誤差項を小さくすることは可能であるが、
同じ計算量であれば4次のルンゲクッタの刻み幅を小さくするほうが精度が良い。


--------------------------------------------------------------------------------