修改rk4方法的C实现

Question

坦率地说，我的问题是，我不确定这是如何工作的。

我需要修改双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 

现在介绍一个dum我的可变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返回v的i组分的t的导数。使用这些名称重新编写的2D系统中，我们得到：

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

这恰恰是功能f一样。

Answer 2

N被定义为常数。这意味着这些循环将要做2次迭代，i = 0和i = 1

f()该函数将返回如果i == 0并通过在多项式的第二元件该多项式的第一个元素的负数，如果i == 1。

我不知道获取谐振子（听起来有点像鹰眼拉法格会说需要重新校准什么，老老实实）的公式，但我会假设，就是这样。