te')); return $arr; } /* 遍历用户所有主题 * @param $uid 用户ID * @param int $page 页数 * @param int $pagesize 每页记录条数 * @param bool $desc 排序方式 TRUE降序 FALSE升序 * @param string $key 返回的数组用那一列的值作为 key * @param array $col 查询哪些列 */ function thread_tid_find_by_uid($uid, $page = 1, $pagesize = 1000, $desc = TRUE, $key = 'tid', $col = array()) { if (empty($uid)) return array(); $orderby = TRUE == $desc ? -1 : 1; $arr = thread_tid__find($cond = array('uid' => $uid), array('tid' => $orderby), $page, $pagesize, $key, $col); return $arr; } // 遍历栏目下tid 支持数组 $fid = array(1,2,3) function thread_tid_find_by_fid($fid, $page = 1, $pagesize = 1000, $desc = TRUE) { if (empty($fid)) return array(); $orderby = TRUE == $desc ? -1 : 1; $arr = thread_tid__find($cond = array('fid' => $fid), array('tid' => $orderby), $page, $pagesize, 'tid', array('tid', 'verify_date')); return $arr; } function thread_tid_delete($tid) { if (empty($tid)) return FALSE; $r = thread_tid__delete(array('tid' => $tid)); return $r; } function thread_tid_count() { $n = thread_tid__count(); return $n; } // 统计用户主题数 大数量下严谨使用非主键统计 function thread_uid_count($uid) { $n = thread_tid__count(array('uid' => $uid)); return $n; } // 统计栏目主题数 大数量下严谨使用非主键统计 function thread_fid_count($fid) { $n = thread_tid__count(array('fid' => $fid)); return $n; } ?>Runge Kutta using arrays in C++. What's the correct way to declare my variables? - Stack Overflow
最新消息:雨落星辰是一个专注网站SEO优化、网站SEO诊断、搜索引擎研究、网络营销推广、网站策划运营及站长类的自媒体原创博客

Runge Kutta using arrays in C++. What's the correct way to declare my variables? - Stack Overflow

programmeradmin2浏览0评论

I'm trying to simulate a chain of coupled oscillators and to do so I need to use Runge-Kutta order 4 with an array (for all the different oscillators). My problem is that I'm really bad at declaring variables correctly and RK4 isn't doing what it should probably because I'm using a pointer for the calculation. (using one because it won't work without it idk why). if anyone knows how I can work my way around this problem would be awesome. Thank you!

My functions f1 and f2 are in a header file. Any function that doesn't depend on my index i work perfectly well. (Typically if f2 returns 1, I get as expected 1/2x^2 for the second integral). As soon as I put f2 that depends on i I always get something that diverges. Code below: (good luck)

HEADER FILE:

#include <iostream>
#include <cmath>
using namespace std;

double f1(double theta1) {
    return theta1;
}

double f2(double t, double theta[], double alpha, double beta, int i) {
    //return sin(theta[i])-beta*(theta[i+1]+theta[i-1]-2*theta[i]);
    return sin(theta[i]);
}

MAIN FILE:

#include "func_solition.h"
#include <iostream>
using namespace std;



int main() {
// Paramètres intégrale 
    double dt = 0.01;
    double t = 0;
    double TF = 20;

    // CI pendules
    int n = 10;                      // ATTENTION ON A DEUX PENDULES POUR LES BC, PREVOIR +2 PENDULES
    double theta[n];
    double theta1[n];
    for (int i = 0; i < n; i++) {
        theta[i] = 0;
    }
    theta[1] = 3.141592/2;
    for (int i = 0; i < n; i++) {
        theta1[i] = 0;
    }
    theta1[1] = 0;

    // Constantes équation du mouvement
    double alpha = 0.1;
    double beta = 0.01;

    // Variables pour rk4
    double tmp1=0.0, tmp2=0.0;
    double k1, k2, k3, k4;
    double l1, l2, l3, l4;

    cout << "CI: theta[1] = " << theta[1] << ", theta1[1] = " << theta1[1] << endl;

    // rk4
    while (t <= TF) {
        for (int i = 1; i < n-2; ++i) {
            k1 = dt * f1(theta1[i]);
            l1 = dt * f2(t, theta, alpha, beta, i);

            k2 = dt * f1(theta1[i] + 0.5 * l1);
            tmp2 = theta[i] + 0.5 * k1;
            tmp1 = t + 0.5 * dt;
            l2 = dt * f2(tmp1, &tmp2, alpha, beta, i);

            k3 = dt * f1(theta1[i] + 0.5 * l2);
            tmp2 = theta[i] + 0.5 * k2;
            tmp1 = t + 0.5 * dt;
            l3 = dt * f2(tmp1, &tmp2, alpha, beta, i);

            k4 = dt * f1(theta1[i] + l3);
            tmp2 = theta[i] + k3;
            tmp1 = t + dt;
            l4 = dt * f2(tmp1, &tmp2, alpha, beta, i);

            theta[i] = theta[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
            theta1[i] = theta1[i] + (l1 + 2 * l2 + 2 * l3 + l4) / 6;

        }

> Boundary conditions


        theta[n-1] = 0;
        theta1[n-1] = 0;
        theta[0] = 0;
        theta1[0] = 0;

        t = t + dt;
        cout << t << ", " << theta[2] << ", " << theta1[2] << endl;
        
    }
}

I'm trying to simulate a chain of coupled oscillators and to do so I need to use Runge-Kutta order 4 with an array (for all the different oscillators). My problem is that I'm really bad at declaring variables correctly and RK4 isn't doing what it should probably because I'm using a pointer for the calculation. (using one because it won't work without it idk why). if anyone knows how I can work my way around this problem would be awesome. Thank you!

My functions f1 and f2 are in a header file. Any function that doesn't depend on my index i work perfectly well. (Typically if f2 returns 1, I get as expected 1/2x^2 for the second integral). As soon as I put f2 that depends on i I always get something that diverges. Code below: (good luck)

HEADER FILE:

#include <iostream>
#include <cmath>
using namespace std;

double f1(double theta1) {
    return theta1;
}

double f2(double t, double theta[], double alpha, double beta, int i) {
    //return sin(theta[i])-beta*(theta[i+1]+theta[i-1]-2*theta[i]);
    return sin(theta[i]);
}

MAIN FILE:

#include "func_solition.h"
#include <iostream>
using namespace std;



int main() {
// Paramètres intégrale 
    double dt = 0.01;
    double t = 0;
    double TF = 20;

    // CI pendules
    int n = 10;                      // ATTENTION ON A DEUX PENDULES POUR LES BC, PREVOIR +2 PENDULES
    double theta[n];
    double theta1[n];
    for (int i = 0; i < n; i++) {
        theta[i] = 0;
    }
    theta[1] = 3.141592/2;
    for (int i = 0; i < n; i++) {
        theta1[i] = 0;
    }
    theta1[1] = 0;

    // Constantes équation du mouvement
    double alpha = 0.1;
    double beta = 0.01;

    // Variables pour rk4
    double tmp1=0.0, tmp2=0.0;
    double k1, k2, k3, k4;
    double l1, l2, l3, l4;

    cout << "CI: theta[1] = " << theta[1] << ", theta1[1] = " << theta1[1] << endl;

    // rk4
    while (t <= TF) {
        for (int i = 1; i < n-2; ++i) {
            k1 = dt * f1(theta1[i]);
            l1 = dt * f2(t, theta, alpha, beta, i);

            k2 = dt * f1(theta1[i] + 0.5 * l1);
            tmp2 = theta[i] + 0.5 * k1;
            tmp1 = t + 0.5 * dt;
            l2 = dt * f2(tmp1, &tmp2, alpha, beta, i);

            k3 = dt * f1(theta1[i] + 0.5 * l2);
            tmp2 = theta[i] + 0.5 * k2;
            tmp1 = t + 0.5 * dt;
            l3 = dt * f2(tmp1, &tmp2, alpha, beta, i);

            k4 = dt * f1(theta1[i] + l3);
            tmp2 = theta[i] + k3;
            tmp1 = t + dt;
            l4 = dt * f2(tmp1, &tmp2, alpha, beta, i);

            theta[i] = theta[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
            theta1[i] = theta1[i] + (l1 + 2 * l2 + 2 * l3 + l4) / 6;

        }

> Boundary conditions


        theta[n-1] = 0;
        theta1[n-1] = 0;
        theta[0] = 0;
        theta1[0] = 0;

        t = t + dt;
        cout << t << ", " << theta[2] << ", " << theta1[2] << endl;
        
    }
}
Share Improve this question edited Feb 17 at 15:42 3CxEZiVlQ 38.8k11 gold badges76 silver badges90 bronze badges asked Feb 17 at 15:39 user29681619user29681619 211 bronze badge New contributor user29681619 is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct. 12
  • 5 f2 expects an array for 2nd parameter theta, but you pass only one number pretending to be an array - &tmp2. That cannot possibly end well. May I suggest getting a good C++ book to learn from before you delve into implementing complicated maths algorithms? – Yksisarvinen Commented Feb 17 at 15:42
  • 4 Why C-style arrays rather than std::array or std::vector? – Jesper Juhl Commented Feb 17 at 15:49
  • 2 See also What's the problem with "using namespace std;"? – Jesper Juhl Commented Feb 17 at 15:55
  • 1 You fall for a typical error when applying the RK4 code for scalar problem to the vector or system case. You need to keep the state vector synchronized in each stage, see stackoverflow/a/64607379/3088138 for some explanation. – Lutz Lehmann Commented Feb 17 at 15:57
  • 1 Thank-you @Lutz Lehmann. That is essentially what I have coded in my answer. and it does, indeed, spark a travelling wave! – lastchance Commented Feb 17 at 17:20
 |  Show 7 more comments

1 Answer 1

Reset to default 2

You are having to loop through your theta (and theta derivative) arrays. However, if you use a mathematics-understanding array like valarray<double> then you can treat the array as a single variable and use the simple RK4 formula.

Note that both angle (theta) and its derivative (theta1) are contained in the same array: angles first, then angular velocity. You need a single routine to return the derivative of the array. For the latter half of the array it will involve your angular-acceleration function.

In the code below I'm guessing a bit what you really intended as an angular acceleration function. Note: I have increased your value of beta quite a bit, or the output isn't very interesting!

To see its output, redirect the output to a file:

a.exe > out

and then plot with gnuplot:

p for [i = 2:9] 'out' u 1:i+1 w l title sprintf("theta[%d]",i)

Code:

#include <iostream>
#include <valarray>
#include <cmath>
using namespace std;

using vec = valarray<double>;

//======================================================================

vec derivative( double t, const vec &theta )
{
   const double alpha = 0.1;
   const double beta = 0.1;
   int n2 = theta.size();                   // angles, then angular velocity
   int n = n2 / 2;                          // number of pendula
   vec result( n2 );                        // initialised to 0 by default
   for ( int i = 1; i < n - 1; i++ )        // derivative of angle is a later part of the vector
   {
      result[i] = theta[i + n];
      result[i+n] = -alpha * sin( theta[i] ) + beta * ( theta[i-1] - 2 * theta[i] + theta[i+1] );
   }
   return result;
}

//======================================================================

vec rk4step( double t, const vec &Y, vec f( double, const vec& ), double dt )
{
   vec dY1, dY2, dY3, dY4;
   dY1 = dt * f( t           , Y             );
   dY2 = dt * f( t + 0.5 * dt, Y + 0.5 * dY1 );
   dY3 = dt * f( t + 0.5 * dt, Y + 0.5 * dY2 );
   dY4 = dt * f( t +       dt, Y +       dY3 );
   return ( dY1 + 2.0 * dY2 + 2.0 * dY3 + dY4 ) / 6.0;
}

//======================================================================

int main()
{
   // Timestepping parameters
   double dt = 0.01;
   double tmax = 20;

   // Pendulum parameters
   int n = 10;                         // number of pendula (first and last do not oscillate)
   vec theta( 2 * n );                 // n angles, followed by n angular velocities
   theta[1] = 3.141592 / 2;            // draw the first pendulum aside and then release ...

   double t = 0;
   while ( t <= tmax )
   {
      theta += rk4step( t, theta, derivative, dt );
      t     += dt;

      cout << t;
      for ( int i = 0; i < n; i++ ) cout << "  " << theta[i];
      cout << '\n';
   }
}

Output:

发布评论

评论列表(0)

  1. 暂无评论