Неявный метод Рунге Кутта для решения жесткой системы ДУ

[Someone]

(Snegovik)

Snegovik в сообщении #1478220 писал(а):

Почему то отступы в коде удваиваются, становятся шире чем в среде разработки.

Они не удваиваются, по прежнему одинарные. Просто шрифты другие, с более широкими пробелами.

[Pphantom]

(Оффтоп)

Snegovik в сообщении #1478220 писал(а):

ps Почему то отступы в коде удваиваются, становятся шире чем в среде разработки.

Разные настройки показа табуляции. в среде - 4 пробела, тут - стандартные 8.

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

[Snegovik]

В коде ошибка, столкнулся в вопросом решения уравнения с одного из форумов:

$\dfrac{dy}{dx} = -1000\cdot y^2$

при $y(0)=10\quad$ на промежутке $x=\left\lbrace0\dots 1\right\rbrace$


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

Вот этот блок кода

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

  error = 0
  ytemp1 = y0
  i=0
       
    do while (abs(g(y, h, x, ytemp2)) > eps)
           dg = (g(y, h, x, ytemp1+dy) - g(y, h, x, ytemp1-dy)) / (dy+dy)
           ytemp2 = y - g(y, h, x, ytemp1) / dg
           ytemp1 = ytemp2
           i = i + 1
              if (i > 10000) then
                 error = 1
                 exit
              end if
    end do

  y2 = ytemp2 
 

Нужно заменить вот этим блоком кода

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

  error = 0
  ytemp1 = y    !Начальное приближение на основе найденного на предыдущем шаге
  !ytemp1 = y0  !При другом начальном приближении
  ytemp2 = 7.33 !Любое значение, далее оно будет определено
  i=0
        
    do while (abs(g(y, h, x, ytemp2)) > eps)
                
      dg = (g(y, h, x, ytemp1+dy) - g(y, h, x, ytemp1-dy)) / (dy+dy) 
      ytemp2 = ytemp1 - g(y, h, x, ytemp1) / dg
      ytemp1 = ytemp2
      i = i + 1
                        
        if (i > 500) then
           error = 1
           exit
        end if
                                                                
    end do
        
  y2 = ytemp2   
 

В коде самой строчки метода вместо $y$ должно было быть $ytemp_1$, как сразу это проглядел непонятно. Плюсом перед циклом сделаны небольшие изменения, назначение ошибки при $500$ итерациях достаточно вместо $10000$. Дополнительно введено определение максимального числа итераций метода Ньютона по любому из обращений к модулю, но это здесь не показано.

ps Теперь вышеприведенное уравнение решается даже с количеством шагов $n=1$ ! правда точность при $n=1$ не достаточно высокая, нужно больше шагов, но не для устойчивости, а для точности.