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

[Lorein_]

Есть закон радиоактивного распада $\frac{du}{dx} = -\lambda x $ С начальным условием $U(0)=1$
Задание: Нужно написать метод РК4 с контролем локальной погрешности.
Сам метод записан.В цикле вызываю его и поточечно строю в chart. Но графиком получается прямая параллельная Ох.Не могу понять где тут ошибка и как правильно нарисовать график.

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Drawing.Drawing2D;
using System.Runtime.InteropServices;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
 
 
namespace RK4_1
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
          
        }
      
 
        double f(double x, double u)
        {
            double L = double.Parse(textBox5.Text);
            return (-L * x);
        }
 
        double rk4(double xn, double yn, double h)
        {
            double k1, k2, k3, k4, yn1;
            k1 = h * f(xn, yn);
            k2 = h * f(xn + 0.5 * h, yn + 0.5 * k1);
            k3 = h * f(xn + 0.5 * h, yn + 0.5 * k2);
            k4 = h * f(xn + h, yn + k3);
            yn1 = yn + (1.0 / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
            return yn1;
        }
 
        private void textBox1_TextChanged(object sender, EventArgs e)
        {
            try
            {
                double x = double.Parse(textBox1.Text);
                errorProvider1.SetError(textBox1, "");
 
            }
            catch
            {
 
                errorProvider1.SetError(textBox1, "Недопустимое значение параметра ");
            }
        }
 
        private void textBox2_TextChanged(object sender, EventArgs e)
        {
            try
            {
                double x = double.Parse(textBox2.Text);
                errorProvider1.SetError(textBox2, "");
 
            }
            catch
            {
 
                errorProvider1.SetError(textBox2, "Недопустимое значение параметра ");
            }
        }
 
        private void textBox3_TextChanged(object sender, EventArgs e)
        {
            try
            {
                double x = double.Parse(textBox3.Text);
                errorProvider1.SetError(textBox3, "");
 
            }
            catch
            {
 
                errorProvider1.SetError(textBox3, "Недопустимое значение параметра a ");
            }
        }
 
        private void textBox4_TextChanged(object sender, EventArgs e)
        {
            try
            {
                double x = double.Parse(textBox4.Text);
                errorProvider1.SetError(textBox4, "");
 
            }
            catch
            {
 
                errorProvider1.SetError(textBox4, "Недопустимое значение параметра b ");
            }
        }
 
        private void textBox7_TextChanged(object sender, EventArgs e)
        {
            try
            {
                double x = double.Parse(textBox7.Text);
                errorProvider1.SetError(textBox7, "");
 
            }
            catch
            {
 
                errorProvider1.SetError(textBox7, "Недопустимое значение параметра h ");
            }
        }
 
        private void textBox5_TextChanged(object sender, EventArgs e)
        {
            try
            {
                double x = double.Parse(textBox5.Text);
                if (x>1 || x< -1) MessageBox.Show("Значение x должно лежать в пределах [0;1]", "Ошибка", MessageBoxButtons.OK, MessageBoxIcon.Error);
            
                errorProvider1.SetError(textBox5, "");
 
            }
            catch
            {
 
                errorProvider1.SetError(textBox5, "Недопустимое значение параметра L ");
            }
        }
 
        private void button2_Click(object sender, EventArgs e)
        {
            double eps= double.Parse(textBox2.Text);
            double x = double.Parse(textBox3.Text);//начало интервала
            double u = double.Parse(textBox6.Text);
            double h = double.Parse(textBox7.Text);//начальный шаг
            double b = double.Parse(textBox4.Text); //конец интервала
            double N = int.Parse(textBox8.Text); //число итераций
            double Vi=0, V_2i=0,S=0,V=0;
 
            chart1.Series[0].Points.Clear();
            dataGridView1.Visible = true;
            dataGridView1.Rows.Clear();
            dataGridView1.Columns.Clear();
            dataGridView1.ColumnCount = 9;
            dataGridView1.Columns[0].HeaderText = "i";
            dataGridView1.Columns[1].HeaderText = "h";
            dataGridView1.Columns[2].HeaderText = "x";
            dataGridView1.Columns[3].HeaderText = "f(x,u)";
            dataGridView1.Columns[4].HeaderText = "V";
            dataGridView1.Columns[5].HeaderText = "Vi";
            dataGridView1.Columns[6].HeaderText = "V_2i";
            dataGridView1.Columns[7].HeaderText = "S";
            for (int i = 0; i < N; i++)
            {
                V = rk4(x,u,h);
                Vi = rk4(x, u, h/2);
                V_2i = rk4(x+h/2, u+h/2, h / 2);
                S = (V_2i - Vi) / 15;
                if ((Math.Abs(S) <= -eps) || (Math.Abs(S) >= eps))
                {
                    h /= 2;
                }
                if (Math.Abs(S) < eps / Math.Pow(2,5))
                {
                    h *= 2;
                }
                if ((Math.Abs(S) > eps / Math.Pow(2, 5)) && (Math.Abs(S) < eps))
                {
                    x += h;
                   
                }
                chart1.Series[0].Points.AddXY(x, rk4(x, u, h));
                dataGridView1.Rows.Add();
                dataGridView1.Rows[i].Cells[0].Value= i;
                dataGridView1.Rows[i].Cells[1].Value = h;
                dataGridView1.Rows[i].Cells[2].Value = x;
                dataGridView1.Rows[i].Cells[3].Value = f(x, u);
                dataGridView1.Rows[i].Cells[4].Value = V;
                dataGridView1.Rows[i].Cells[5].Value = Vi;//  значение c половинным шагом
                dataGridView1.Rows[i].Cells[6].Value = V_2i; //половинный шаг в точке x+h/2
                dataGridView1.Rows[i].Cells[7].Value = S;
 
                
 
               x += h;
 
            }
            
 
 
 
        }
 
 
 
     
    }
}

[Pphantom]

Выкиньте для начала эту рисовальную шелуху и проверьте, что численное решение правильное.

P.S. И если вы привели задание дословно (т.е. в нем не сказано явно, что нужно рисовать графики в той же программе) - рисовальную шелуху нужно выкинуть не для начала, а насовсем. Ввод данных - из текстового файла (если он вообще зачем-то нужен), вывод результатов - в другой текстовый файл.

[guryev]

Lorein_ в сообщении #1426572 писал(а):

закон радиоактивного распада $\frac{du}{dx} = -\lambda x $

Не похожа эта штука на радиоактивный распад. Может быть, так: $\frac{du}{dx} = -\lambda u$?

[Lorein_]

Да,конечно $\frac{dU}{dx} = -\lambda u$. Очепятка вышла)

[Lorein_]

Pphantom в сообщении #1426581 писал(а):

Выкиньте для начала эту рисовальную шелуху и проверьте, что численное решение правильное.

P.S. И если вы привели задание дословно (т.е. в нем не сказано явно, что нужно рисовать графики в той же программе) - рисовальную шелуху нужно выкинуть не для начала, а насовсем. Ввод данных - из текстового файла (если он вообще зачем-то нужен), вывод результатов - в другой текстовый файл.

Нет,задание приведено не дословно.У меня задача стоит сделать программу с интерфейсом решающую оду методом рунге-кутта 4-го порядка.
Метод проверила,и видимо в нём и есть ошибка.Но я не могу понять где именно.Формулы записаны верно

[Pphantom]

Lorein_ в сообщении #1426804 писал(а):

У меня задача стоит сделать программу с интерфейсом решающую оду методом рунге-кутта 4-го порядка.

Печально.

Lorein_ в сообщении #1426804 писал(а):

Метод проверила,и видимо в нём и есть ошибка.Но я не могу понять где именно.Формулы записаны верно

Т.е. численные результаты не соответствуют ожидаемым?

[Lorein_]

Pphantom в сообщении #1426808 писал(а):

Lorein_ в сообщении #1426804 писал(а):

У меня задача стоит сделать программу с интерфейсом решающую оду методом рунге-кутта 4-го порядка.

Печально.

Lorein_ в сообщении #1426804 писал(а):

Метод проверила,и видимо в нём и есть ошибка.Но я не могу понять где именно.Формулы записаны верно

Т.е. численные результаты не соответствуют ожидаемым?

Да.Проверяю на ду $ \frac{dy}{dx}= 2x-y+x^2, y(0)=0 $. Его точное решение $x^2$. Однако в таблице в этой колонке числа далёкие от правды получаются.

[Pphantom]

Ну в собственно методе Рунге-Кутты все и в самом деле верно. А вот с правилом Рунге вы творите что-то странное (что хотелось получить вот в этом участке кода:

Используется синтаксис C#

V = rk4(x,u,h);
Vi = rk4(x, u, h/2);
V_2i = rk4(x+h/2, u+h/2, h / 2);
S = (V_2i - Vi) / 15;

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

То, что у вас все входные данные нормально парсятся (и потом выходные нормально выводятся), проверяли? Попробуйте просто задать параметры для счета в коде (без ввода из всяких формочек) и результаты в файл вывести.

[EUgeneUS]

А где в цикле $u = rk4(x, u, h)$ ?

Или я чего-то не понимаю...

-- 19.11.2019, 21:21 --

                if ((Math.Abs(S) > eps / Math.Pow(2, 5)) && (Math.Abs(S) < eps))
                {
                    x += h;
                   
                }

Это выглядит лишним.
Выше в зависимости от $S$ корректируем размер шага. А тут, хоба! зачем-то на шаг сдвинулись.

[Lorein_]

Для оценки локальной погрешности мне необходимо вычислить $S = \frac{\hat{u_{n+1}}-u_{n+1}}{2^p -1}$
где $\hat{u_{n+1}}$ - значение в точке $(x_{n+h/2},u_{n+h/2})$ c шагом $h/2$,$ u_{n+1}$ -значение в точке $(x_n,u_n)$ c шагом h . Соответственно порядок метода 4-й.В знаменателе получаем 15. В коде V - численное решение,
Vi - соответствует $u_{n+1}$, V_2i - $\hat{u_{n+1}}$.
Получается,если заданное $|S|< \varepsilon$ то шаг в два раза меньше. Если $|S|<\frac{ \varepsilon}{2^p -1}$ - шаг удваиваем. Если $\frac{ \varepsilon}{2^p -1}<|S|<\varepsilon$ - продолжаем с тем же шагом
Если задать параметры для счёта в коде,то ничего не меняется.Считает также неправильно.
Я сомневаюсь в цикле.Правильно ли вообще это записано? Численно решаю диффур и увеличиваю x на h.

for (int i = 0; i < n; i++)
            {
                V = rk4(x, u, h);
                chart1.Series[0].Points.AddXY(x, rk4(x, u, h));
                dataGridView1.Rows.Add();
                dataGridView1.Rows[i].Cells[0].Value = i;
                dataGridView1.Rows[i].Cells[1].Value = x;
                dataGridView1.Rows[i].Cells[2].Value = h;
                dataGridView1.Rows[i].Cells[3].Value = f(x, u);
                dataGridView1.Rows[i].Cells[4].Value = V;
                dataGridView1.Rows[i].Cells[5].Value = Vi;//  значение c половинным шагом
                dataGridView1.Rows[i].Cells[6].Value = V_2i; //половинный шаг в точке x+h/2
                dataGridView1.Rows[i].Cells[7].Value = S;
               x += h;

            }

Цитата:

А где в цикле $u = u + rk4(x, u, h)$ ?

Или я чего-то не понимаю...

Если добавить,то например для x=0 считает V = 1,38179331667133.

[EUgeneUS]

Lorein_ в сообщении #1426822 писал(а):

Если добавить

В конец цикла надо добавлять, там где Вы $x$ увеличиваете.
Иначе в каждой итерации цикла вы подаете на вход $rk4( x, u, h)$ начальное значение $u$, а не текущее.

И чуть другое надо добавить - я исправил свой пост выше.

[Lorein_]

Цитата:

Выше в зависимости от $S$ корректируем размер шага. А тут, хоба! зачем-то на шаг сдвинулись.

И правда!Чушь какая-то)
И всё равно значения получаются не такими.

[Lorein_]

Цитата:

В конец цикла надо добавлять, там где Вы $x$ увеличиваете.

Так как точное решение оду $ \frac{dy}{dx}= 2x-y+x^2, y(0)=0 $, $x^2$ то и численное решение должно быть похоже на квадрат x.
Например на промежутке [0;2] с шагом h=0.2. Однако в программе расчёты получаются совершенно другими. Не могу понять в чем же может быть дело.

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;

namespace RungeKutta
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }

        double f(double x, double u)
        {
                      return 2*x-u+x*x;
        }

         double rk4(double x,double u,double h) {
            double k1, k2, k3, k4;
            //double b = double.Parse(textBox6.Text); //конец интервала
            double b = 2.2;
            double n=(b-x)/h;
          for(int i = 0; i < n; i++)
            {
                k1 = h*f(x,u);
                k2 = h * f(x + h / 2, u + (k1)/ 2);
                k3 = h * f(x + h / 2, u + (k2)/ 2);
                k4 = h *f(x + h, u + k3);
                u=u+ (k1 + 2 * k2 + 2 * k3 + k4)*(h / 6);
                x += h;  
            }
            return u;
        }

        private void button1_Click(object sender, EventArgs e)
        {
         
            double x = 0;
            double u = 0;
            double h = 0.2;
            double b = 2.2;

            double n = (b - x) / h;
          
            double Vi = 0, V_2i = 0, S = 0, V = 0;

            chart1.Series[0].Points.Clear();
            dataGridView1.Visible = true;
            dataGridView1.Rows.Clear();
            dataGridView1.Columns.Clear();
            dataGridView1.ColumnCount = 9;
            dataGridView1.Columns[0].HeaderText = "i";
            dataGridView1.Columns[1].HeaderText = "x";
            dataGridView1.Columns[2].HeaderText = "h";
            dataGridView1.Columns[3].HeaderText = "f(x,u)";
            dataGridView1.Columns[4].HeaderText = "V";
            dataGridView1.Columns[5].HeaderText = "Vi";
            dataGridView1.Columns[6].HeaderText = "V_2i";
            dataGridView1.Columns[7].HeaderText = "S";
            for (int i = 0; i < n; i++)
            {
                V = rk4(x, u, h);
                chart1.Series[0].Points.AddXY(x, rk4(x, u, h));
                dataGridView1.Rows.Add();
                dataGridView1.Rows[i].Cells[0].Value = i;
                dataGridView1.Rows[i].Cells[1].Value = x;
                dataGridView1.Rows[i].Cells[2].Value = h;
                dataGridView1.Rows[i].Cells[3].Value = f(x, u);
                dataGridView1.Rows[i].Cells[4].Value = V;
                dataGridView1.Rows[i].Cells[5].Value = Vi;
                dataGridView1.Rows[i].Cells[6].Value = V_2i; 
                dataGridView1.Rows[i].Cells[7].Value = S;
               x += h;
                u = rk4(x, u, h);
                
            }

        }

[EUgeneUS]

У Вас функция $rk4$ что делает?
Считает следующую итерацию? А зачем там цикл появился?
Оставьте её как была в начале.

(Оффтоп)

Генетическое программирование какое-то.

[Lorein_]

Спасибо за помощь!)