Метод Рунге-Кутты 4 порядка на 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] + " " + "|"); } } }

Ответ

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