K
K
Kari JR2015-09-15 10:44:24
Programming
Kari JR, 2015-09-15 10:44:24

Euler task, what is the function F responsible for in this code?

Gentlemen, a lot of questions arose about the Euler method and the implementation of the numerical integration of a system of differential equations. It was possible with grief in half to dig up the source code for solving this problem, but unfortunately not everything is clear in it, where I figured out I put comments. But there are still a lot of ambiguities for me. For example, what the function f0, f1, etc. is responsible for. It would be great if you could point out errors, if any, to a more correct solution, and comment on the places that raise questions. Thank you!

//---------------------------------------------------------------------------

#pragma hdrstop
#include <stdio.h>
#include <math.h>
//---------------------------------------------------------------------------

#pragma argsused
double f0(double t);
double f1(double t, double y[]);
double f2(double t, double y[]);
double f3(double t, double y[]);
void eiler( double t, double y[], int n, double tn, double tk, int m, double delt);
int main(int argc, char* argv[])
{
        int n; //порядок системы дифференциальных уравнений
        double delt; //шаг интегрирования
        double tn; //начальное значение интервала интегрирования;//
        double tk; //конечное значение интервала интегрирования
        int m ; //шаг интегрирования
        double y[3];
        double ys[3];
        double t;
        int j;
        int i;
        printf("Input n="); //Порядок интегрирования
        scanf("%d" , &n);
        if(n  > 3)
        {
                return-1;
        }
        printf("Input delt=");
        scanf("%lf" , &delt);
        printf("Input tn=");
        scanf("%lf" , &tn);
        printf("Input tk=");
        scanf("%lf" , &tk);
        m = (tk-tn)/delt;
   
           //Считываем начальные данные
        for ( j= 0; j <n; j++)
        {
              printf("Input y^(%d",j);
              printf(") =");
              scanf("%lf" , &y[j]);

        }
           // решение

                t = tn;
        for ( j= 0; j <=m; j++)
        {
              eiler(t, y,  n, tn,tk, m,  delt);
              t = t + delt;
              printf("x= %lf", t);
              for (i = 0; i <n ; i++)
              {
                printf("   y^%d", i);
                printf("= %lf", y[i]);
              }
              printf("\n");

        }

        scanf("%s" , &n);
        return 0;
}
//---------------------------------------------------------------------------
double f0(double t)
{
   return(5+t*t);
}
//---------------------------------------------------------------------------
double f1(double t,double  y[])
{
   return(5+t*t-y[0]);
}
//---------------------------------------------------------------------------
double f2(double t, double  y[])
{
   return(5+t*t-3*y[1]-y[0]);
}
//---------------------------------------------------------------------------
double f3(double t, double  y[])
{
   return(5+t*t-2*y[2]-3*y[1]-y[0]);
}


 void eiler( double t, double y[], int n, double tn, double tk, int m, double delt)
{

      switch (n )
      {
        case 3:
        {
        //Порядок 3
                y[2] = y[2] + delt * f3(t, y);
                y[1] = y[1] + delt * y[2];
                y[0] = y[0] + delt * y[1];
                break;
        }
        case 2:
        {
        //Порядок 2
                y[1] = y[1] + delt * f2(t, y);
                y[0] = y[0] + delt * y[1];
                break;
        }
        case 1:
        {
        //Порядок 1
               y[0] = y[0] + delt * f1(t, y);
                break;
        }
        default:
        {  
    //Порядок непонятен выводим 0
                y[0] = f0(t);
                break;
        }
}
}

Answer the question

In order to leave comments, you need to log in

1 answer(s)
M
Mercury13, 2015-09-15
@KariJR

Качество кода — на уровне студенческой лабораторной. Причём очень хреновой.
Программа решает четыре системы.
y = f0(t)
y' = f1(t, y)
y'' = f2(t, y, y')
y''' = f3(t, y, y'')
1. Если уж говорить о методе Эйлера — в программе ошибка.
Тут нужно старое значение y[2], а не новое.
2. Функцию eiler стоило бы обозвать eulerStep, заодно убрав лишние параметры.
3. Автор не знает такой концепции, как процедурный тип aka callback, поэтому работа с правой частью у него вышла вот через такую задницу.
4. В промышленном коде метод решения ОДУ я бы вынес в отдельную процедуру, с callback’ами на правую часть и на потребителя результата. Причём внешнее (доступное пользователю библиотеки) именование должно быть такое, чтобы её пользователь мог ею воспользоваться, даже если он не читал статью, по которой библиотеку писали. Так что смотрите, какие имена являются стандартизированными, понятными по любому учебнику, а какие стоит прояснить. Внутреннее именование — это уж думайте сами, научный код обычно сложен, и его постоянно приходится сличать со статьёй. Поэтому имена как в статье — для научного кода не WTF.
5. Не слишком удачно разделена ответственность между функцией правой части и функцией шага по Эйлеру. Это привело к ограничению — программа может решать только ОДУ y{n} = f(t, y, y', ..., y{n−1}).

Didn't find what you were looking for?

Ask your question

Ask a Question

731 491 924 answers to any question