14-QT5 GSL求解微分方程并结合qcustomplot绘图

最新推荐文章于 2024-02-29 20:12:40 发布

落叶_小唱

最新推荐文章于 2024-02-29 20:12:40 发布

阅读量1.2k

点赞数 1

分类专栏： Qt Math 文章标签： GSL qcustomplot

本文链接：https://blog.csdn.net/ouening/article/details/100613236

版权

Math 同时被 2 个专栏收录

37 篇文章 3 订阅

订阅专栏

22 篇文章 3 订阅

订阅专栏

1. 所求微分方程

$2y′−y=4sin(3t),\quad y(0)=y0,\quad t \in [0,8]$

2. GSL相关函数

定义ODE系统：

   gsl_odeiv2_system
   
   // 所求函数
   int (* function) (double t, const double y[], double dydt[], void * params)
   // 所求函数的雅克比形式函数
   int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params)
   // 微分方程阶数（纬度）
   size_t dimension
   // 其他参数
   void * params

以前面定义求解的微分方程为例，写作如下代码（需要一阶导数前面的系数需要化为1）：

 #include "math.h"
    int func (double t, const double y[], double f[],
          void *params)
    {
      (void)(t); /* avoid unused parameter warning */
      double mu = *(double *)params;
      f[0] = 2*sin(3*t)+y[0]/2;
      return GSL_SUCCESS;
    }

可以计算出雅克比矩阵函数：

 int jac (double t, const double y[], double *dfdy,
         double dfdt[], void *params)
 {
      (void)(t); /* avoid unused parameter warning */
   double mu = *(double *)params;
      gsl_matrix_view dfdy_mat
        = gsl_matrix_view_array (dfdy, 1, 1);
      gsl_matrix * m = &dfdy_mat.matrix;
      gsl_matrix_set (m, 0, 0, 0.5);
    
      dfdt[0] = 0.0;
        
      return GSL_SUCCESS;
    }

定义求解方法（driver）

GSL中求解微分方程的方法有很多种，有常规步长的方法，有自适应步长控制的方法（Adaptive Step-size Control）以及进化方法（evolution），GSL支持的gsl_odeiv2_step_type有如下(参考GSL官方帮助手册第357页）：

  - gsl_odeiv2_step_rk2
  - gsl_odeiv2_step_rk4
  - gsl_odeiv2_step_rkf45
    
  - gsl_odeiv2_step_rkck
    
  - gsl_odeiv2_step_rk8pd
    
  - gsl_odeiv2_step_rk1imp
    
  - gsl_odeiv2_step_rk2imp
    
  - gsl_odeiv2_step_rk4imp

  - gsl_odeiv2_step_bsimp

  - gsl_odeiv2_step_msadams

  - gsl_odeiv2_step_msbdf

应用求解

int gsl_odeiv2_driver_apply(gsl_odeiv2_driver * d, double * t, const double t1, double y[])
- 参数1 *d：微分方程求解方法，gsl_odeiv2_driver类型
- 参数2 *t：初始时刻
- 参数3 t1：微分方程需要计算的时刻，该值是改变的
- 参数4 y[]：含初始时刻t的初值，如果是二阶微分方程，其形式类似y[]={1,0}

3. qcustomplot绘图

自己编写一个函数，实现微分方程的求解以及绘图操作（省去了相关头文件及参数定义）

void MainWindow::setupQuadraticDemo(QCustomPlot *customPlot)
{
    QString demoName = "GSL ODE Demo";
    int len = 300;  // data length
    // generate some data:
    QVector<double> x(len), y(len); // initialize with entries 0..100

    // ---------------GLS ODE-------------------//
    double mu = NULL;
    gsl_odeiv2_system sys = {func, jac, 1, &mu};

    gsl_odeiv2_driver * d =
            gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
                                           1e-6, 1e-6, 0.0);
    double t = 0.0, b = 10.0;
    double a = t;
    double y0[1] = { -24.0/37.0 };       // initial value at t is 1.0

    for (int i = 0; i <= len; i++)
    {
        double h = (b-a)/(double)len;
        double ti = a + i * h;
        int status = gsl_odeiv2_driver_apply (d, &t, ti, y0);

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

    gsl_odeiv2_driver_free (d);
    //====================END GSL ODE=====================

    // create graph and assign data to it:
    customPlot->addGraph();
    customPlot->graph(0)->setData(x, y);
    // give the axes some labels:
    customPlot->xAxis->setLabel("x");
    customPlot->yAxis->setLabel("y");
    // set axes ranges, so we see all data:
    customPlot->xAxis->rescale(true);
    customPlot->yAxis->rescale(true);
}

4. mainwindow.h完整代码

#ifndef MAINWINDOW_H
#define MAINWINDOW_H

#include <QMainWindow>
#include "qcustomplot.h"	// 导入qcustomplot

namespace Ui {
class MainWindow;
}

class MainWindow : public QMainWindow
{
    Q_OBJECT

public:
    explicit MainWindow(QWidget *parent = nullptr);
    ~MainWindow();

    void setupQuadraticDemo(QCustomPlot *customPlot);

private:
    Ui::MainWindow *ui;
};

#endif // MAINWINDOW_H

5. mainwindow.cpp完整代码：

#include "mainwindow.h"
#include "ui_mainwindow.h"
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_sf_bessel.h>
#include "math.h"

// df/dt = 2sin(3t)+y/2
int func (double t, const double y[], double f[],
          void *params)
{
    (void)(t); /* avoid unused parameter warning */
    double mu = *(double *)params;
    f[0] = 2*sin(3*t)+y[0]/2;
    //    f[0] = t-1+1/t-2*y[0];
    return GSL_SUCCESS;
}

int jac (double t, const double y[], double *dfdy,
         double dfdt[], void *params)
{
    (void)(t); /* avoid unused parameter warning */
    double mu = *(double *)params;
    gsl_matrix_view dfdy_mat
            = gsl_matrix_view_array (dfdy, 1, 1);
    gsl_matrix * m = &dfdy_mat.matrix;
    gsl_matrix_set (m, 0, 0, 0.5);
    //     gsl_matrix_set (m, 0, 0, -2);

    //    dfdt[0] = 0.0;
    dfdt[1] = -1.0;

    return GSL_SUCCESS;
}


MainWindow::MainWindow(QWidget *parent) :
    QMainWindow(parent),
    ui(new Ui::MainWindow)
{
    ui->setupUi(this);

    setupQuadraticDemo(ui->customPlot);
    setWindowTitle("QCustomPlot: GSL ODE Demo");
    statusBar()->clearMessage();
    //    currentDemoIndex = demoIndex;
    ui->customPlot->replot();
}

MainWindow::~MainWindow()
{
    delete ui;
}

void MainWindow::setupQuadraticDemo(QCustomPlot *customPlot)
{
    int len = 300;  // data length
    // generate some data:
    QVector<double> x(len), y(len); // initialize with entries 0..100

    // ---------------GLS ODE-------------------//
    double mu = NULL;
    gsl_odeiv2_system sys = {func, jac, 1, &mu};

    gsl_odeiv2_driver * d =
            gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
                                           1e-6, 1e-6, 0.0);
    double t = 0.0, b = 10.0;
    double a = t;
    double y0[1] = { -24.0/37.0 };       // initial value at t is 1.0

    for (int i = 0; i <= len; i++)
    {
        double h = (b-a)/(double)len;
        double ti = a + i * h;
        int status = gsl_odeiv2_driver_apply (d, &t, ti, y0);

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

    gsl_odeiv2_driver_free (d);
    //====================END GSL ODE=====================

    // create graph and assign data to it:
    customPlot->addGraph();
    customPlot->graph(0)->setData(x, y);
    // give the axes some labels:
    customPlot->xAxis->setLabel("x");
    customPlot->yAxis->setLabel("y");
    // set axes ranges, so we see all data:
    customPlot->xAxis->rescale(true);
    customPlot->yAxis->rescale(true);
    //    customPlot->xAxis->setRange(-1, t1);
    //    customPlot->yAxis->setRange(-3, 3);
    //    customPlot->yAxis->scaleRange(1);
}

6. 注意事项

导入qcustomplot之后pro文件需要支持printsupport:QT += widgets printsupport

pro文件需要导入GSL头文件路径以及库路径，以Linux为例：

INCLUDEPATH += /usr/local/include
LIBS += -L/usr/local/lib -lgsl -lgslcblas -lm

mainwindow.ui界面拖入一个widget控件，将其提升为QCustomPlot类，将控件名改为customPlot
程序组织结构
实例结果显示：
-

落叶_小唱

关注

1
点赞
踩
13

收藏

觉得还不错? 一键收藏
2
评论
14-QT5 GSL求解微分方程并结合qcustomplot绘图

1. 所求微分方程2y′−y=4sin(3t),y(0)=y0,t∈[0,8]2y′−y=4sin(3t),\quad y(0)=y0,\quad t \in [0,8]2y′−y=4sin(3t),y(0)=y0,t∈[0,8]2. GSL相关函数定义ODE系统： gsl_odeiv2_system // 所求函数 int (* function) (doub...
复制链接

扫一扫

专栏目录