Модификация C реализации метода rk4

Question

Моя проблема, честно говоря, в том, что я не уверен, как это работает.

Мне нужно изменить функцию double f (), чтобы решить произвольное дифференциальное уравнение d2θ / dt2 = −ω2sinθ, но так какЯ не уверен, что делать дальше.

Я понимаю саму функцию rk4 runge4 ();Что я не понимаю, так это то, как функция f () возвращает правильные значения для гармонического осциллятора.

Может кто-нибудь хотя бы объяснить логику функции f ()?

Оригиналкод ниже.

/* 
************************************************************************
*  rk4.c: 4th order Runge-Kutta solution for harmonic oscillator       *
*                      *
* From: "A SURVEY OF COMPUTATIONAL PHYSICS" 
   by RH Landau, MJ Paez, and CC BORDEIANU 
   Copyright Princeton University Press, Princeton, 2008.
   Electronic Materials copyright: R Landau, Oregon State Univ, 2008;
   MJ Paez, Univ Antioquia, 2008; & CC BORDEIANU, Univ Bucharest, 2008
   Support by National Science Foundation                              
*
************************************************************************
*/

#include <stdio.h>

#define N 2                                   /* number of equations */
#define dist 0.1                              /* stepsize */
#define MIN 0.0                               /* minimum x */
#define MAX 10.0                              /* maximum x */

void runge4(double x, double y[], double step);
double f(double x, double y[], int i);

int main()
{

   double x, y[N];
   int j;

   FILE *output;                              /* save data in rk4.dat */
   output = fopen("rk4.dat","w");

   y[0] = 1.0;                                /* initial position  */
   y[1] = 0.0;                                /* initial velocity  */

   fprintf(output, "%f\t%f\n", x, y[0]);

   for(x = MIN; x <= MAX ; x += dist)
   {
      runge4(x, y, dist);
      fprintf(output, "%f\t%f\n", x, y[0]);   /* position vs. time */
   }
   printf("data stored in rk4.dat\n");
   fclose(output);
}
/*-----------------------end of main program--------------------------*/

/* Runge-Kutta subroutine */
void runge4(double x, double y[], double step)
{
   double h=step/2.0,                         /* the midpoint */
          t1[N], t2[N], t3[N],                /* temporary storage */
          k1[N], k2[N], k3[N],k4[N];          /* for Runge-Kutta  */
   int i;

   for (i=0; i<N; i++) t1[i] = y[i]+0.5*(k1[i]=step*f(x, y, i));
   for (i=0; i<N; i++) t2[i] = y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
   for (i=0; i<N; i++) t3[i] = y[i]+    (k3[i]=step*f(x+h, t2, i));
   for (i=0; i<N; i++) k4[i] =                 step*f(x + step, t3, i);

   for (i=0; i<N; i++) y[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
/*--------------------------------------------------------------------*/

/* definition of equations - this is the harmonic oscillator */
double  f(double x, double y[], int i)
{
   if (i == 0) return(y[1]);               /* RHS of first equation */
   if (i == 1) return(-y[0]);              /* RHS of second equation */
}

Answer 1

Начните с закона Гука:

F = -kx

Объедините это со вторым законом Ньютона, чтобы получить дифференциальное уравнение для линейного гармонического осциллятора:

ma = F = -kx
mx'' = -kx
x'' = -k/m x

Произвольно выбрал наши единицы так, чтобы k/m == 1, и уравнение становится просто:

x'' = -x

Теперь введите фиктивную переменную y = x' и запишите это дифференциальное уравнение второго порядка в виде двумерной системы первого порядка:

x' = y
y' = -x

Функция f в вашем коде кодирует именно эту систему; Я собираюсь изменить имена переменных для ясности:

double  f(double t, double v[], int i)
{
   if (i == 0) return(v[1]);
   if (i == 1) return(-v[0]);
}

v - это вектор [x,y] из вышеприведенной двумерной системы. Для i, t и v функция f возвращает производную по t i-го компонента v. Переписав 2d систему под этими именами, мы получим:

dv[0]/dt =  v[1]
dv[1]/dt = -v[0]

Именно это и делает функция f.

Answer 2

N определяется как константа 2. Это означает, что эти циклы будут выполнять 2 итерации, i = 0 и i = 1

Функция f() вернет второй элемент переданного полинома, если i == 0, и отрицательный результат первого элемента этого полинома, если i == 1.

Я не знаю формулу для получения гармонического осциллятора (звучит так, как сказал бы Джорди Лафорж, нуждается в повторной калибровке или что-то в этом роде), но я бы предположил, что это так.