#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
/* Рунге-Кутта для X'=f(t,X) X - вектор.
Параметры: t, x - массив X, xs - массив X',
n - размерность, f- ф-ция-вектор - возвращает рез-ты в массив Y
f(t,x,Y)
mem - память 5*n*sizeof(double) или NULL...
*/
int RKn(int n, double*t, double* x, double h,
void (*f)(int, double, double*, double*), double * mem)
{
double *k1,*k2,*k3,*X,*Y;
int j;
if (mem) k1 = mem;
else if ((k1
= (double*)malloc(5*sizeof(double)*n
))==0) return 0;
Y = (X = (k3 = (k2 = k1+n)+n)+n)+n;
f(n,*t,X,Y);
for(j=0; j<n; j++) X[j] = x[j] + (k1[j]=h*Y[j])/2.;
f(n,*t+h/2.,X,Y);
for(j=0; j<n; j++) X[j] = x[j] + (k2[j]=h*Y[j])/2.;
f(n,*t+h/2.,X,Y);
for(j=0; j<n; j++) X[j] = x[j] + (k3[j]=h*Y[j]);
f(n,*t+=h,X,Y);
for(j=0; j<n; j++) x[j]+=(k1[j]+2.*(k2[j]+k3[j])+h*Y[j])/6.;
return 1;
};
void f(int n, double t, double*x, double*Y)
{
Y
[0] = -52*x
[0]-100*x
[1]+exp(-t
);}
int main()
{
double t = 0.0;
double x[2] = {1.0, 0.0};
double h = 0.01;
printf("%4.2lf %12.7lf %12.7lf\n",t
,x
[0],x
[1]);
while(t < 2.0)
{
if (!RKn(2,&t,x,h,f,0))
{
}
printf("%4.2lf %12.7lf %12.7lf\n",t
,x
[0],x
[1]); }
return 0;
}
I2luY2x1ZGUgPHN0ZGxpYi5oPgojaW5jbHVkZSA8c3RkaW8uaD4KI2luY2x1ZGUgPG1hdGguaD4KI2luY2x1ZGUgPHN0cmluZy5oPgoKLyog0KDRg9C90LPQtS3QmtGD0YLRgtCwINC00LvRjyBYJz1mKHQsWCkgWCAtINCy0LXQutGC0L7RgC4KICAg0J/QsNGA0LDQvNC10YLRgNGLOiB0LCB4IC0g0LzQsNGB0YHQuNCyIFgsIHhzIC0g0LzQsNGB0YHQuNCyIFgnLAogICBuIC0g0YDQsNC30LzQtdGA0L3QvtGB0YLRjCwgZi0g0YQt0YbQuNGPLdCy0LXQutGC0L7RgCAtINCy0L7Qt9Cy0YDQsNGJ0LDQtdGCINGA0LXQty3RgtGLINCyINC80LDRgdGB0LjQsiBZCiAgIGYodCx4LFkpCiAgIG1lbSAtINC/0LDQvNGP0YLRjCA1Km4qc2l6ZW9mKGRvdWJsZSkg0LjQu9C4IE5VTEwuLi4KKi8KCmludCBSS24oaW50IG4sIGRvdWJsZSp0LCBkb3VibGUqIHgsIGRvdWJsZSBoLAogICAgICAgIHZvaWQgKCpmKShpbnQsIGRvdWJsZSwgZG91YmxlKiwgZG91YmxlKiksIGRvdWJsZSAqIG1lbSkKewogICAgZG91YmxlICprMSwqazIsKmszLCpYLCpZOwogICAgaW50IGo7CiAgICBpZiAobWVtKSBrMSA9IG1lbTsKICAgIGVsc2UgaWYgKChrMSA9IChkb3VibGUqKW1hbGxvYyg1KnNpemVvZihkb3VibGUpKm4pKT09MCkgcmV0dXJuIDA7CgogICAgWSA9IChYID0gKGszID0gKGsyID0gazErbikrbikrbikrbjsKICAgIG1lbWNweShYLHgsc2l6ZW9mKGRvdWJsZSkqbik7CgogICAgZihuLCp0LFgsWSk7CiAgICBmb3Ioaj0wOyBqPG47IGorKykgWFtqXSA9IHhbal0gKyAoazFbal09aCpZW2pdKS8yLjsKICAgIGYobiwqdCtoLzIuLFgsWSk7CiAgICBmb3Ioaj0wOyBqPG47IGorKykgWFtqXSA9IHhbal0gKyAoazJbal09aCpZW2pdKS8yLjsKICAgIGYobiwqdCtoLzIuLFgsWSk7CiAgICBmb3Ioaj0wOyBqPG47IGorKykgWFtqXSA9IHhbal0gKyAoazNbal09aCpZW2pdKTsKICAgIGYobiwqdCs9aCxYLFkpOwogICAgZm9yKGo9MDsgajxuOyBqKyspIHhbal0rPShrMVtqXSsyLiooazJbal0razNbal0pK2gqWVtqXSkvNi47CiAgICBmcmVlKGsxKTsKCiAgICByZXR1cm4gMTsKfTsKCgp2b2lkIGYoaW50IG4sIGRvdWJsZSB0LCBkb3VibGUqeCwgZG91YmxlKlkpCnsKICAgIFlbMF0gPSAtNTIqeFswXS0xMDAqeFsxXStleHAoLXQpOwogICAgWVsxXSA9IHhbMF0gKyBzaW4odCk7Cn0KCmludCBtYWluKCkKewogICAgZG91YmxlIHQgPSAwLjA7CiAgICBkb3VibGUgeFsyXSA9IHsxLjAsIDAuMH07CiAgICBkb3VibGUgaCA9IDAuMDE7CgogICAgcHJpbnRmKCIlNC4ybGYgICAgJTEyLjdsZiAgICAlMTIuN2xmXG4iLHQseFswXSx4WzFdKTsKCiAgICB3aGlsZSh0IDwgMi4wKQogICAgewogICAgICAgIGlmICghUktuKDIsJnQseCxoLGYsMCkpCiAgICAgICAgewogICAgICAgICAgICBwcmludGYoIkVycm9yIVxuIik7CiAgICAgICAgfQogICAgICAgIHByaW50ZigiJTQuMmxmICAgICUxMi43bGYgICAgJTEyLjdsZlxuIix0LHhbMF0seFsxXSk7CiAgICB9CiAgICByZXR1cm4gMDsKfQo=