Рунге-Кутта 2-го порядка дифференциальных уравнений


Я мог бы как-то улучшить эту (не численно, только c-мудрый)? Мне нужно, чтобы выделить массив, так как я ее вернуть, а не просто использовать его внутри функции?

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Approximates a solution to a differential equation on the form: 
   y'(t) + ay(t) = x(t)
   y(0) = b 
*/
double* runge_kutta_2nd_order(double stepSize, double a, double b, double (*x) (double), double upto)
{
    int resultSize = ((int) (upto / stepSize)) + 2;
    double yt = b;
    double time;
    double k1,k2,ystar1,ystar2;
    int index = 1;

    //double *results = (double*) malloc(resultSize * (sizeof(double)));
    double *results = malloc(resultSize * sizeof(*results));
    if(results == NULL)
    {
        printf("\nCould not allocate memory. Exiting program.");
        exit(0);
    }

    results[0] = b;

    for(time = 0; time < upto; time += stepSize) //<=
    {   
        k1 = x(time) - a * yt;
        ystar1 = yt + stepSize * k1;
        k2 = x(time + stepSize) - a * ystar1;
        ystar2 = yt + (k1 + k2) / 2 * stepSize;
        yt = ystar2;
        results[index] = ystar2;
        index++;
    }
    return results;
}

void free_results(double **r)
{
    free(*r);
    *r = NULL;
}


double insignal(double t)
{
    return exp(t/2)*(sin(5*t) - 10*cos(5*t));
}

int main(void)
{
    int i;
    double *res = runge_kutta_2nd_order(0.01,-1,0,&insignal,10);


    printf("\nRunge Kutta 2nd order approximation of the differential equation:");
    printf("\ny'(t) - y(t) = e^(t/2) * (sin(5t) - 10cos(5t))");
    printf("\ny(0) = 0");
    printf("\n0 <= t <= 10");

    for(i=0; i<1001; i+=100){
        printf("\ni = %lf => y = ", 0.01*i);
        printf("%lf", res[i]);
    }
    printf("\n");

    free_results(&res);


    return 0;
}


337
2
задан 4 декабря 2011 в 08:12 Источник Поделиться
Комментарии
1 ответ

Это больше обычного, чтобы взять указатель на уже выделенный массив, вместо того чтобы выделить массив внутри функции. Это позволяет вызывающей функции, чтобы выбрать лучший вид распределения.

При этом вы, скорее всего, также возьмите количество индексов для заполнения в качестве входных данных, расчета размера шаге от этого:

/* resultsize must be >= 2 */
void runge_kutta_2nd_order(double results[], size_t resultsize, double a, double b, double (*x) (double), double upto)
{
const double stepSize = upto / (resultsize - 1);
double yt = b;
size_t index;

results[0] = b;

for (index = 1; index < resultsize; index++)
{
const double time = index * stepSize;
const double k1 = x(time) - a * yt;
const double ystar1 = yt + stepSize * k1;
const double k2 = x(time + stepSize) - a * ystar1;
const double ystar2 = yt + (k1 + k2) / 2 * stepSize;
yt = ystar2;
results[index] = ystar2;
}
}

Так, например, в этом случае можно использовать автоматическое распределение в функции main():

int main(void)
{
int i;
double res[1001];

runge_kutta_2nd_order(res, 1001, -1, 0, &insignal, 10);

printf("\nRunge Kutta 2nd order approximation of the differential equation:");
printf("\ny'(t) - y(t) = e^(t/2) * (sin(5t) - 10cos(5t))");
printf("\ny(0) = 0");
printf("\n0 <= t <= 10");

for(i=0; i<1001; i+=100){
printf("\ni = %f => y = ", 0.01*i);
printf("%f", res[i]);
}
printf("\n");

return 0;
}

3
ответ дан 5 декабря 2011 в 12:12 Источник Поделиться