Необходимо Решить задачу Коши для обыкновенного дифференциального уравнения, в компиляторе численного решения выдает NAN

Рейтинг: 0Ответов: 1Опубликовано: 04.07.2023

Решить задачу Коши для обыкновенного дифференциального уравнения y'= f(x,y) на промежутке [a b, ] методом Рунге – Кутты четвертого порядка точности: y'= (2*x -5)/x^2 *y+5 y(a)=b , a=2, b=4 Точное решение задачи: y(x)=x^2

Код программы

Не могу понять, почему делится на 0

#include <math.h>

#include <conio.h>

#include <stdio.h>

#define n 100

float f(float, float);

float ya(float);

int main (void)

{

 float a, b, h, x, y[n];

 float m1,m2,m3,m4; // дополнительные переменные

 int k;

 
 a=2; b=4; h=(b-a)/(n-1); 

 y[2]=4;  // начальное условие

 for (k=0; k<n-1; k++) // метод Рунге – Кутты

 {

  x=k*h;

  m1=h*f(x,y[k]);

  m2=h*f(x+h/2,y[k]+m1/2);

  m3=h*f(x+h/2,y[k]+m2/2);

  m4=h*f(x+h,y[k]+m3);

  y[k+1]=y[k]+(m1+2*m2+2*m3+m4)/6;

 }

 
 // сравнение с аналитическим решением

 for (k=0; k<n; k++) 

 {

  x=k*h;

  printf("numerical = %5.2f \t\t analytical = %5.2f\n",y[k],ya(x));

 }

 getch();

}


float f(float x, float y) // правая часть

{

 return (2*x-5)/pow(x, 2)*y + 5;  

}


float ya(float x) // аналитическое решение

{

 return pow(x, 2); 

}
c# nan численные-методы
Источник: Stack Overflow на русском

Ответы

▲ 2Принят

Я бы делал так, без всяких массивов, и с правильной инициализацией переменных. Кстати, у вас не С++, а чистый С.

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

#define n 100

double f(double, double);

double ya(double);

int main(void)
{
    double a = 2, b = 4, h = (b-a)/n, x = a, y = b;

    printf("%7.3lf   %10.6lf    %10.6lf\n",x,y,ya(x));
    while(x < b) // метод Рунге-Кутты
    {

        double m1 = h*f(x,y);
        double m2 = h*f(x+h/2,y+m1/2);
        double m3 = h*f(x+h/2,y+m2/2);
        double m4 = h*f(x+h,y+m3);
        x += h;
        y += (m1 + 2*(m2+m3) + m4)/6;

        printf("%7.3lf   %10.6lf    %10.6lf\n",x,y,ya(x));
    }
}


double f(double x, double y) // правая часть

{
    return (2*x-5)/(x*x)*y+5;
}


double ya(double x) // аналитическое решение
{
    return x*x;
}