常微分方程边值问题求解

原创
2023/12/05 11:01
阅读数 29

最近需要求解一个常微分方程边值问题如下:

\(\frac{d^2u}{dx^2}-ab(u^4-c)=0\)

其边条件如下:

\(u(0)=T_0, \qquad u'(L)=-a[u^4(L)-c]\)

式中a, b, c, T0和L均为常数。

常规的求解方式,是把高阶微分方程化为一阶常微分方程组,然后求解。但这通常用于对初值问题的求解,对于边值问题,则需要稍加改进。一般是采用初值问题叠加一个寻根方法,来找到合适的边值。考虑到该方程关于x是反射对称的,因此把方程化简为:

\(\frac{du}{dx}=v\)

\(\frac{dv}{dx}=ab(u^4-c)\)

其边条件也化简为:

\(v(0)=au^4(0)-ac\)

\(u(L)=T_0\)

如此,则可定义如下微分方程初值问题:

\(\frac{du}{dx}=v\)

\(\frac{dv}{dx}=ab(u^4-c)\)

\(v(0)=au^4(0)-ac\)

\(u(0)=T\)

求解该微分方程组,可得到u(x, T),其中T为试探参数。接下来考虑一维寻根问题的解:

\(u(x, T)|_{x=L} - T_0 = 0\)

如此即可得到最开始边值问题的解u(x)。

 

下面给出代码,计算用到了GSL库,其windows版本可在如下地址获得:

https://github.com/BrianGladman/gsl

其他相关版本可参考如下链接处的说明。

https://www.gnu.org/software/gsl/extras/native_win_builds.html

 

test.cpp

#define WIN32
#define _CRT_SECURE_NO_WARNINGS
#define GSL_DLL

#include <stdio.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_errno.h>


/**********************************************************/
/* ODE test system definitions                            */
/**********************************************************/

struct parameters {
    double a;
    double b;
    double c;
    double L;
    double T0;
};

int func(double t, const double y[], double f[], void* params)
{
    struct parameters* p = (struct parameters*)params;

    f[0] = y[1];
    f[1] = p->a * p->b * y[0] * y[0] * y[0] * y[0] - p->a * p->b * p->c;

    return GSL_SUCCESS;
}

int jac(double t, const double y[], double* dfdy, double dfdt[], void* params)
{
    struct parameters* p = (struct parameters*)params;

    dfdy[0] = 0.0;
    dfdy[1] = 1.0;

    dfdy[2] = 4 * p->a * p->b * y[0] * y[0] * y[0];
    dfdy[3] = 0.0;

    dfdt[0] = 0.0;
    dfdt[1] = 0.0;

    return GSL_SUCCESS;
}

double diff_sys(double T, void* params)
{
    struct parameters* p = (struct parameters*)params;

    gsl_odeiv2_system sys = {
      func,
      jac,
      2,
      p
    };

    gsl_odeiv2_driver* d =
        gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd,
            1e-6, 1e-6, 0.0);

    double t0 = 0.0, tL = p->L;
    double y[2];
    y[0] = T;
    y[1] = p->a * y[0] * y[0] * y[0] * y[0] - p->a * p->c;

    for (int i = 1; i <= 100; i++)
    {
        double ti = i * tL / 100.0;
        int status = gsl_odeiv2_driver_apply(d, &t0, ti, y);

        if (status != GSL_SUCCESS)
        {
            printf("error, return value=%d\n", status);
            break;
        }

        printf("%.5e %.5e %.5e\n", t0, y[0], y[1]);
    }

    gsl_odeiv2_driver_free(d);
    return y[0] - p->T0;
}

/**********************************************************/
/* Main function                                          */
/**********************************************************/

int main(void) 
{
    struct parameters p;
    p.a = 0.85e-10;
    p.b = 4.0;
    p.c = 16.0;
    p.T0 = 100;
    p.L = 1.0;


    int status;
    int iter = 0, max_iter = 100;
    const gsl_root_fsolver_type* T;
    gsl_root_fsolver* s;
    double r = 0, r_expected = 100;
    double x_lo = p.T0 - 100.0, x_hi = p.T0 + 100.0;
    gsl_function F;

    F.function = &diff_sys;
    F.params = &p;

    T = gsl_root_fsolver_brent;
    s = gsl_root_fsolver_alloc(T);
    gsl_root_fsolver_set(s, &F, x_lo, x_hi);

    printf("using %s method\n",
        gsl_root_fsolver_name(s));

    printf("%5s [%9s, %9s] %9s %10s %9s\n",
        "iter", "lower", "upper", "root",
        "err", "err(est)");

    do
    {
        iter++;
        status = gsl_root_fsolver_iterate(s);
        r = gsl_root_fsolver_root(s);
        x_lo = gsl_root_fsolver_x_lower(s);
        x_hi = gsl_root_fsolver_x_upper(s);
        status = gsl_root_test_interval(x_lo, x_hi,
            1e-3, 1e-3);

        if (status == GSL_SUCCESS)
            printf("Converged:\n");

        printf("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
            iter, x_lo, x_hi,
            r, r - r_expected,
            x_hi - x_lo);
    } while (status == GSL_CONTINUE && iter < max_iter);

    gsl_root_fsolver_free(s);

    return status;
}

 

完整的工程可在如下地址下载:

https://download.csdn.net/download/u014559935/88599245

展开阅读全文
加载中
点击引领话题📣 发布并加入讨论🔥
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部