2015-01-13 71 views
2

我正在研究使用GSL的非線性微分方程。事情是我對C的東西很新。我只是將GNU網站上的示例調整到了我現在感興趣的等式中。如何諮詢帶許多參數和諧波函數的GSL ODE

這是公式:

d2x/dt2 + r*dx/dy + cos(x) + v*cos(2*x+0.4) E1*sin(wt) + E2*sin(2*w*t+a) = 0 

什麼我堅持是我不知道如何在代碼多個參數堵塞。此外,我不知道如何在此代碼中使用餘弦或正弦函數。

我試圖找出這個問題,通過在谷歌一路搜索。我找不到任何有助於我的事情。

#include <stdio.h> 
#include <gsl/gsl_errno.h> 
#include <math.h> 
#include <gsl/gsl_math.h> 
#include <gsl/gsl_matrix.h> 
#include <gsl/gsl_odeiv2.h> 

int func (double t, const double x[], double y[], void *params) 
{ 
    double r = *(double *)params; 
    double v = *(double *)params; 
    double w = *(double *)params; 
    double E1 = *(double *)params; 
    double E2 = *(double *)params; 
    double a = *(double *)params; 
    y[0] = x[1]; 
    y[1] = -r*x[1] - cos(x[0]) - v*cos(2*x[0]+0.4) - E1*sin(w*t) - E2*sin(2*w*t+a); 
    return GSL_SUCCESS; 
} 

int jac (double t, const double x[], double *dydx, double dydt[], void *params) 
{ 
    double r = *(double *)params; 
    double v = *(double *)params; 
    double w = *(double *)params; 
    double E1 = *(double *)params; 
    double E2 = *(double *)params; 
    double a = *(double *)params; 
    gsl_matrix_view dydx_mat = gsl_matrix_view_array (dydx, 2, 2); 
    gsl_matrix * m = &dydx_mat.matrix; 
    gsl_matrix_set (m, 0, 0, 0.0); 
    gsl_matrix_set (m, 0, 1, 1.0); 
    gsl_matrix_set (m, 1, 0, sin(x[0]) + 2*v*sin(2*x[0]+0.4)); 
    gsl_matrix_set (m, 1, 1, -r); 
    dydt[0] = 0.0; 
    dydt[1] = 0.0; 
    return GSL_SUCCESS; 
} 

int main (void) 
{ 
    double r = 0.0; 
    double v = 0.0; 
    double w = 2.4; 
    double E1 = -2.3; 
    double E2 = 0; 
    double a = 0.7; 
    gsl_odeiv2_system sys = {func, jac, 2, &r, &v, &w, &E1, &E2, &a}; 

    gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_x_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0); 
    int i; 
    double t = 0.0, t1 = 10000; 
    double x[2] = {0.0, 0.0}; 

    for (i = 1 ; i<=10000; i++) 
     { 
      double ti = i*t1/10000; 
      int status = gsl_odeiv2_driver_apply (d, &t, ti, x); 

      if (status != GSL_SUCCESS) 
       { 
        printf("error, return value%d\n", status); 
        break; 
       } 
      printf("%.5e %.5e %.5e\n", t, x[0], x[1]); 
     } 

    gsl_odeiv2_driver_free (d); 
    return 0; 
} 
+1

你確定你輸入的公式正確嗎?你有一個'dx/dy'的術語。 – Rufflewind

回答

4

params參數是一個指針(地址/存儲位置),以一些任意數據結構。在來自GSL文檔的示例中,它們的公式僅包含一個參數,這意味着只傳遞double-精確數字的地址即可。

但是,對於您的問題,您需要訪問6個不同的參數。您無法訪問具有相同地址的每個參數!

/* this doesn't work! */ 
double r = *(double *)params; 
double v = *(double *)params; 
double w = *(double *)params; 
double E1 = *(double *)params; 
double E2 = *(double *)params; 
double a = *(double *)params; 

由於所有的地址都相同,所以您指的是相同的號碼。要解決這個問題,您可以:將所有參數存儲在長度爲6的數組中,或將它們存儲在預定義的數據結構中。後一種方法更具可讀性,所以我將證明這一點。

首先定義數據類型指定什麼參數,你將存儲:

struct param_type { 
    double r; 
    double v; 
    double w; 
    double E1; 
    double E2; 
    double a; 
}; 

現在,創建這種類型的main功能的結構和存儲參數的實際值:

struct param_type my_params = {r, v, w, E1, E2, a}; 

當定義系統時,您存儲一個指針struct param_type

gsl_odeiv2_system sys = {func, jac, 2, &my_params}; 

要使用參數中funcjac,只需從一個普通的指針(void *)投了params參數指針爲特定的數據類型(struct param_type *):

struct param_type *my_params_pointer = params; 

(請注意,在C++這必須可以與明確投寫的)最後,您可以通過訪問參數:

double r = my_params_pointer->r; 
double v = my_params_pointer->v; 
double w = my_params_pointer->w; 
double E1 = my_params_pointer->E1; 
double E2 = my_params_pointer->E2; 
double a = my_params_pointer->a; 

箭頭->這裏代替的點.,因爲my_params_pointer指針並且在使用前需要爲dereferenced