Метод Рунге-Кутты 4 порядка на java

#java #методы #численные_методы

Я только начала изучать java, и многие элементарные вещи для меня являются непонятными.

Прощу вас подсказать мне, как решать систему взаимозависимых уравнений методом Рунге-Кутты!

Полное задание:

Для решения полученной задачи Коши для системы первого порядка вида:

y'= f(t,y), y(0)=y0

использовать метод Рунге_Кутты 4-го порядка точности:

k1 = f(tn , yn)
 k2 = f(tn + h/2 , yn + hk1 / 2)
 k3 = f(tn + h/2 , yn + hk2 / 2)
 k4 = f(tn + h , yn + h*k3)

yn+1=yn+h*(k1 + 2*k2 + 2*k3 + k4) / 6    

На отрезке [0,5] с точным решением 

y1=cos(x)/(1+e2x)1/2,

y2=sin(x)/(1+e2x)1/2.

Для проверки правильности работы программы решить тестовую задачу из двух уравнений:

y1' = –y2 + y1(y12 + y22 – 1),

y2' = y1 + y2(y12 + y22 – 1),

На отрезке [0,5] с точным решением 

y1=cos(x)/(1+e2x)1/2,

y2=sin(x)/(1+e2x)1/2.

Я пытаюсь реализовать только метод!
Сложность заключается в том, что производные у1 и у2 зависят друг от друга.
И правильно ли  я сделала, что производную в точке вычисляю по у1 и у2

public static double y1=Math.cos(x)/Math.pow(1 + Math.pow(Math.E, 2 * x),0.5);
public static double y2=Math.sin(x)/Math.pow(1+Math.pow(Math.E,2*x),0.5);

//   dy1/dx
public static double derviY1(double x,double y10,double y20){

    return -y20+y10*(Math.pow(y10,2)+Math.pow(y20,2)-1);
}

//   dy2/dx
public static  double derviY2(double x ,double y1,double y2){
    return y1 + y2*(Math.pow(y1,2) + Math.pow(y2,2) - 1);
}


И почему-то вычисляя y20[i+1] и y10[i+1] они у меня остаются неизменными, хотя я
в программе изменяю данные, которые в них входят.

for (int i = 0; i < n - 1; i++) {
    x = i * h;

    k1 = h * derviY1(x, y10[i], y20[i]);
    m1 = h * derviY2(x, y10[i], y20[i]);

    k2 = h * derviY1(x + h / 2, y10[i] + k1 / 2, y20[i] + k1 / 2);
    m2 = h * derviY2(x + h / 2, y10[i] + m1 / 2, y20[i] + m1 / 2);

    k3 = h * derviY1(x + h / 2, y10[i] + k2 / 2, y20[i] + k2 / 2);
    m3 = h * derviY2(x + h / 2, y10[i] + m2 / 2, y20[i] + m2 / 2);

    k4 = h * derviY1(x + h, y10[i] + k3, y20[i] + k3);
    m4 = h * derviY2(x + h, y10[i] + m3, y20[i] + m3);

    y10[i + 1] = y10[i] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
    y20[i + 1] = y20[i] + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6;

    System.out.println("| " + x + " |" + " " + y10[i] + " " + "|" + " " + y20[i]
+ " " + "|");


Полный код программы:

import java.*;
import java.lang.Math.*;
import static java.lang.System.out;

public class Test {

    static double x;
    //start
    private static int a = 0;
    // stop
    private static int b = 5;

    public static double y1 = Math.cos(x) / Math.pow(1 + Math.pow(Math.E, 2 * x), 0.5);
    public static double y2 = Math.sin(x) / Math.pow(1 + Math.pow(Math.E, 2 * x), 0.5);

    //   dy1/dx
    public static double derviY1(double x, double y10, double y20) {

        return -y20 + y10 * (Math.pow(y10, 2) + Math.pow(y20, 2) - 1);
    }

    //   dy2/dx
    public static double derviY2(double x, double y1, double y2) {
        return y1 + y2 * (Math.pow(y1, 2) + Math.pow(y2, 2) - 1);
    }


    public static void main(String[] args) {
        int n = 10;
        double h = (b - a) / n;
        double k1, k2, k3, k4, m1, m2, m3, m4;
        double[] y10 = new double[n]; //array of values y1
        double[] y20 = new double[n]; //array of values y2
        y10[0] = 1 / Math.sqrt(2);
        y20[0] = 0;

        // Computation by 4th order Runge-Kutta
        //update x
        for (int i = 0; i < n - 1; i++) {
            x = i * h;

            k1 = h * derviY1(x, y10[i], y20[i]);
            m1 = h * derviY2(x, y10[i], y20[i]);

            k2 = h * derviY1(x + h / 2, y10[i] + k1 / 2, y20[i] + k1 / 2);
            m2 = h * derviY2(x + h / 2, y10[i] + m1 / 2, y20[i] + m1 / 2);

            k3 = h * derviY1(x + h / 2, y10[i] + k2 / 2, y20[i] + k2 / 2);
            m3 = h * derviY2(x + h / 2, y10[i] + m2 / 2, y20[i] + m2 / 2);

            k4 = h * derviY1(x + h, y10[i] + k3, y20[i] + k3);
            m4 = h * derviY2(x + h, y10[i] + m3, y20[i] + m3);

            y10[i + 1] = y10[i] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
            y20[i + 1] = y20[i] + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6;

            System.out.println("| " + x + " |" + " " + y10[i] + " " + "|" + " " +
y20[i] + " " + "|");
        }
    }
}

    

Ответы

Ответ 1

У вас h = 0. Почему? Ведь вы ясно написали, что h = (5-0)/10. А проблема вот в чём.
Ваше выражение (5-0)/10 сначала преобразуется в int, а потом уже в double. То есть
происходит всё как-то так (int)((5-0)/10). После этой операции мантисса отбрасывается
и получается полноценный 0.
Достаточно привести к double одну из переменных, как всё выражение будет обработано,
как для переменных с плавающей точкой. К примеру так:

double h = ((double)b - a) / n;


В этих преобразованиях есть более интересные моменты, о которых хотелось бы рассказать,
то это уже другая история.