# 常微分方程边值问题求解

2023/12/05 11:01

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

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

$\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)|_{x=L} - T_0 = 0$

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;
}