$\left\{\begin{matrix} y^{''}(x)+p(x)y^{'}(x)+q(x)y(x)=f(x), \space\space a<x<b, \space\space\space\space (1)\\ y(a)=\alpha, \space\space\space\space y(b)=\beta \end{matrix}\right.$

$\left\{\begin{matrix} y^{''}(x)+p(x)y^{'}(x)+q(x)y(x)=f(x), \space\space, a<x<b, \space\space\space\space (2)\\ y^{'}(a)+\lambda y(a)=\alpha, \space\space y^{'}(b)+\mu y(b)=\beta \end{matrix}\right.$

即形如式（1）和（2）的问题计算。其中 $p(x),q(x),f(x)$ 为已知函数， $y(x)$ 为待求函数， $\alpha,\beta,\lambda,\mu$ 为常数。

公式（1）中的边值条件称为狄利克雷边界条件（Dirichlet），公式（2）中的边值条件称为混合边界条件，当 $\lambda=\mu=0$ 时，称为诺依曼边界条件（Neumann）。

二、主要方法及算例实现

1、打靶法

基本思想：将公式（1）转化为与之等价的初值问题，利用假设的待定系数来确定这种转化。

1.1 狄利克雷边值问题

1.1.1 计算思路

先将公式（1）转化为等价的初值问题：

$\left\{\begin{matrix} y^{''}+p(x)y^{'}+q(x)y=f(x),\space\space\space\space (3)\\ y(a)=\alpha, \space\space y^{'}(a)=\gamma \end{matrix}\right.$

其中 $\gamma$ 为常数，待定。为了确定 $\gamma$ ，需要对公式（3）进行分解，该初值问题的解可以通过两个齐次方程和一个非齐次方程来确定，如下：

$\left\{\begin{matrix} y^{''}+p(x)y^{'}+q(x)y=0, \space\space a<x<b,\\ y(a)=1, \space\space y^{'}(a)=0 \end{matrix}\right.$

$\left\{\begin{matrix} y^{''}+p(x)y^{'}+q(x)y=0, \space\space a<x<b,\\ y(a)=0, \space\space y^{'}(a)=1 \end{matrix}\right.$

$\left\{\begin{matrix} y^{''}+p(x)y^{'}+q(x)y=f(x),\space\space a<x<b,\\ y(a)=0, \space\space y^{'}(a)=0 \end{matrix}\right.$

假设上述三式的解分别为 $y_{1}(x),y_{2}(x),y_{3}(x)$ ，根据线性微分方程解的结构可知公式（3）的解为 $y=\alpha y_{1}+\gamma y_{2}+y_{3}$ 。由于公式（1）与公式（3）等价，故公式（3）的解应该满足公式（1）的边界条件， $y(a)=\alpha$ 自然满足，还需满足 $\beta=y(b)=\alpha y_{1}(b)+\gamma y_{2}(b)+y_{3}(b)$ 。于是得到：

$\gamma = \frac{\beta-\alpha y_{1}(b)-y_{3}(b)}{y_{2}(b)}, \space\space\space\space (4)$

于是，求解的基本流程为：

①将边值问题分解为3个等价的初值问题，得到3个数值解，初值问题的解法可以任选（通常采用四阶龙格-库塔法，确保计算精度），且假设这3个数值解为 $y^{1},y^{2},y^{3}$ ，分别对应精确解为 $y_{1},y_{2},y_{3}$ ，精度为四阶。

②用龙格-库塔法求解得到的是 $y^{j}(j=1,2,3)$ 在各节点 $x_{i}(i=0,1,\cdot \cdot \cdot ,m)$ 处的值，记作 $y^{j}_{i}$ ，于是待定系数 $\gamma =\frac{\beta - \alpha y^{1}_{m}-y^{3}_{m}}{y^{2}_{m}}$ 。

③利用 $y_{i}=\alpha y^{1}_{i}+\gamma y^{2}_{i}+y^{3}_{i},i=0,1,\cdot \cdot \cdot ,m$ 得到公式（3），也是公式（1）的精确解 $y(x)$ 在节点 $x_{i},i=0,1,\cdot \cdot \cdot ,m$ 处的数值解。

1.1.2 算例实现

求解两点边值问题：

$\left\{\begin{matrix} y^{''}+y=-1, \space\space 0<x<1,\\ y(0)=y(1)=0 \end{matrix}\right.$

步长h=0.1。已知精确解为： $y=cosx+\frac{1-cos1}{sin1}sinx-1$ 。

代码如下：


#include <cmath>
#include <stdlib.h>
#include <stdio.h>
 
int N;
double h;
 
int main(int argc, char* argv[])
{
        int i;
        double a,b,alpha,beta,gamma;
        double *x,*y,*y1,*y2,*y3;
        double p(double x);
        double q(double x);
        double f(double x);
        double *rhsf(int k, double x, double y1, double y2);
        double *RKhigher(int k, double *x, double c, double d);
        double exact(double x);
 
        a=0.0;
        b=1.0;
        N=10;
        h=(b-a)/N;
        x=(double *)malloc(sizeof(double)*(N+1));
        for(i=0;i<=N;i++)
            x[i]=a+i*h;
 
        alpha=0.0;beta=0.0;
 
        y1=RKhigher(0,x,1,0);
        y2=RKhigher(0,x,0,1);
        y3=RKhigher(1,x,0,0);
 
        y=(double *)malloc(sizeof(double)*(N+1));
        gamma=(beta-alpha*y1[N]-y3[N])/y2[N];
        for(i=1;i<=N;i++)
        {
                y[i]=alpha*y1[i]+gamma*y2[i]+y3[i];
                printf("x[%d]=%.2f, y=%f, exact=%f, error=%.4e.\n",i,x[i],y[i],exact(x[i]),fabs(exact(x[i])-y[i]));
        }
 
        return 0;
}
 
double p(double x)
{
        return 0.0;
}
 
double q(double x)
{
        return 1.0;
}
 
double f(double x)
{
        return -1.0;
}
 
double *rhsf(int k, double x, double y1, double y2)
{
        double *z;
        z=(double *)malloc(sizeof(double)*2);
        z[0]=y2;
        z[1]=k*f(x)-p(x)*y2-q(x)*y1;
        return z;
}
 
double *RKhigher(int k, double *x, double c, double d)
{
        int i;
        double *y,y1,y2,*k1,*k2,*k3,*k4;
        y=(double*)malloc(sizeof(double)*(N+1));
        y[0]=c;
        y1=c;
        y2=d;
        for(i=0;i<N;i++)
        {
                k1=rhsf(k,x[i],y1,y2);
                k2=rhsf(k,x[i]+0.5*h,y1+0.5*h*k1[0],y2+0.5*h*k1[1]);
                k3=rhsf(k,x[i]+0.5*h,y1+0.5*h*k2[0],y2+0.5*h*k2[1]);
                k4=rhsf(k,x[i+1],y1+h*k3[0],y2+h*k3[1]);
                y1=y1+h*(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6.0;
                y2=y2+h*(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6.0;
                y[i+1]=y1;
                free(k1);free(k2);free(k3);free(k4);
        }
        return y;
}
 
double exact(double x)
{
        return cos(x)+(1-cos(1.0))*sin(x)/sin(1.0)-1;
}

N=10时，运行结果如下:


x[1]=0.10, y=0.049543, exact=0.049543, error=8.9716e-08.
x[2]=0.20, y=0.088600, exact=0.088600, error=1.6175e-07.
x[3]=0.30, y=0.116780, exact=0.116780, error=2.1458e-07.
x[4]=0.40, y=0.133801, exact=0.133801, error=2.4707e-07.
x[5]=0.50, y=0.139494, exact=0.139494, error=2.5845e-07.
x[6]=0.60, y=0.133801, exact=0.133801, error=2.4836e-07.
x[7]=0.70, y=0.116780, exact=0.116780, error=2.1683e-07.
x[8]=0.80, y=0.088600, exact=0.088600, error=1.6430e-07.
x[9]=0.90, y=0.049543, exact=0.049543, error=9.1615e-08.
x[10]=1.00, y=-0.000000, exact=0.000000, error=5.5511e-17.

1.2 混合边值问题

1.2.1 计算思路

我们采取同样的思路，将其转化为等价的初值问题，如：

$\left\{\begin{matrix} y^{''}+p(x)y^{'}+q(x)y=f(x),\space\space a<x<b,\\ y(a)=s,\space\space\space\space\space,y^{'}(a)=t \end{matrix}\right.$

其中s，t为常数，待定。 $y=sy_{1}+ty_{2}+y_{3}$ 是上式的解，按照等价条件可以确定：

$s=\frac{\alpha(y^{'}_{2}(b)+ \mu y_{2}(b))-\beta +y^{'}_{3}(b)+\mu y_{3}(b)}{\lambda (y^{'}_{2}(b)+\mu y_{2}(b))-(y^{'}_{1}(b)+\mu y_{1}(b))}$

$t=\frac{\lambda \beta-\alpha (y^{'}_{1}(b)+\mu y_{1}(b))-\lambda (y^{'}_{3}(b)+\mu y_{3}(b))}{\lambda (y^{'}_{2}(b)+\mu y_{2}(b))-(y^{'}_{1}(b)+\mu y_{1}(b))}$

于是，求解的基本流程为：

①降阶化为一阶方程组并分别数值求解常微分方程的初值问题，得到的数值解分别为 $y^{1},z^{1}/y^{2},z^{2}/y^{3},z^{3}$ 。由于二阶初值问题通过引入 $z=y^{'}$ 降阶后就可以将 $y^{'}$ 数值求解，从而得到对应的各导数值的近似值 $z^{1},z^{2},z^{3}$ 。

②数值解 $y^{1},z^{1}/y^{2},z^{2}/y^{3},z^{3}$ 都是离散的，第一步得到的是 $y^{j},z^{j}(j=1,2,3)$ 在各点 $x_{i}(i=0,1,\cdot \cdot \cdot ,m)$ 处的值，分别记作 $y^{j}_{i},z^{j}_{i}$ ，于是待定常数s，t为：

$s=\frac{\alpha (z^{2}_{m}+\mu y^{2}_{m})-\beta+z^{3}_{m}+\mu y^{3}_{m}}{\lambda(z^{2}_{m}+\mu y^{2}_{m})-(z^{1}_{m}+\mu y^{1}_{m})}$

$t=\frac{\lambda \beta-\alpha(z^{1}_{m}+\mu y^{1}_{m})-\lambda(z^{3}_{m}+\mu y^{3}_{m})}{\lambda(z^{2}_{m}+\mu y^{2}_{m})-(z^{1}_{m}+\mu y^{1}_{m})}$

③确定s，t后，再利用 $y_{i}=sy^{1}_{i}+ty^{2}_{m}+y^{3},i=0,1,\cdot \cdot \cdot ,m$ 得到 $y(x)$ 在节点 $x_{i},i=0,1,\cdot \cdot \cdot ,m$ 处的数值解。

1.2.2 算例实现

求解两点边值问题：

$\left\{\begin{matrix} y^{''}+xy^{'}-xy=(x+2)e^{x}cosx, \space\space 0<x<\pi \\ y^{'}(0)-2y(0)=1, \space\space y^{'}(\pi)+y(\pi)=e^{-\pi} \end{matrix}\right.$

步长h=π/10。已知精确解为： $y=e^{x}sinx$ 。

代码如下：


#include <cmath>
#include <stdlib.h>
#include <stdio.h>
#define PI 3.14159265359
 
int N;
double h;
 
int main(int argc, char* argv[])
{
        int i;
        double a,b,alpha,beta,lammda,mu,s,t,eta1,eta2,eta3,temp;
        double *x,*y,*y1,*y2,*y3;
        double p(double x);
        double q(double x);
        double f(double x);
        double *rhsf(int k, double x, double y1, double y2);
        double *RKhigher(int k, double *x, double c, double d);
        double exact(double x);
 
        a=0.0;
        b=PI;
        N=10;
        h=(b-a)/N;
        x=(double *)malloc(sizeof(double)*(N+1));
        for(i=0;i<=N;i++)
            x[i]=a+i*h;
 
        lammda=-2.0;mu=1.0;alpha=1.0;beta=-exp(PI);
        y1=RKhigher(0,x,1,0);
        y2=RKhigher(0,x,0,1);
        y3=RKhigher(1,x,0,0);
        eta1=y1[2*N+1]+mu*y1[N];
        eta2=y2[2*N+1]+mu*y2[N];
        eta3=y3[2*N+1]+mu*y3[N];
        temp=lammda*eta2-eta1;
        s=(alpha*eta2-beta+eta3)/temp;
        t=(lammda*beta-alpha*eta1-lammda*eta3)/temp;
        y=(double*)malloc(sizeof(double)*(N+1));
        for(i=0;i<=N;i++)
        {
                y[i]=s*y1[i]+t*y2[i]+y3[i];
                printf("x[%d]=%.2f, y=%f. exact=%f, error=%.4e.\n",i,x[i],y[i],exact(x[i]),fabs(exact(x[i])-y[i]));
        }
 
        return 0;
}
 
 
double p(double x)
{
        return x;
}
double q(double x)
{
        return -x;
}
double f(double x)
{
        return (x+2)*exp(x)*cos(x);
}
double *rhsf(int k, double x, double y1, double y2)
{
        double *z;
        z=(double*)malloc(sizeof(double)*2);
        z[0]=y2;
        z[1]=k*f(x)-p(x)*y2-q(x)*y1;
        return z;
}
double *RKhigher(int k, double *x, double c, double d)
{
        int i;
        double *y,*z,*w,y1,y2,*k1,*k2,*k3,*k4;
        y=(double*)malloc(sizeof(double)*(N+1));
        z=(double*)malloc(sizeof(double)*(N+1));
 
        y[0]=c;
        z[0]=d;
        y1=c;
        y2=d;
        for(i=0;i<N;i++)
        {
                k1=rhsf(k,x[i],y1,y2);
                k2=rhsf(k,x[i]+0.5*h,y1+0.5*h*k1[0],y2+0.5*h*k1[1]);
                k3=rhsf(k,x[i]+0.5*h,y1+0.5*h*k2[0],y2+0.5*h*k2[1]);
                k4=rhsf(k,x[i+1],y1+h*k3[0],y2+h*k3[1]);
                y1=y1+h*(k1[0]+k2[0]*2+k3[0]*2+k4[0])/6.0;
                y2=y2+h*(k1[1]+k2[1]*2+k3[1]*2+k4[1])/6.0;
                y[i+1]=y1;
                z[i+1]=y2;
                free(k1);free(k2);free(k3);free(k4);
        }
 
        w=(double*)malloc(sizeof(double)*(2*N+2));
        for(i=0;i<=N;i++)
        {
                w[i]=y[i];
                w[N+1+i]=z[i];
        }
        free(y);free(z);
        return w;
}
double exact(double x)
{
        return exp(x)*sin(x);
}

N=10时，运行结果如下：


x[0]=0.00, y=0.001062. exact=0.000000, error=1.0622e-03.
x[1]=0.31, y=0.424819. exact=0.423078, error=1.7411e-03.
x[2]=0.63, y=1.104218. exact=1.101778, error=2.4400e-03.
x[3]=0.94, y=2.079446. exact=2.076207, error=3.2394e-03.
x[4]=1.26, y=3.345909. exact=3.341619, error=4.2907e-03.
x[5]=1.57, y=4.816274. exact=4.810477, error=5.7967e-03.
x[6]=1.88, y=6.271685. exact=6.263717, error=7.9679e-03.
x[7]=2.20, y=7.305872. exact=7.294929, error=1.0943e-02.
x[8]=2.51, y=7.271028. exact=7.256376, error=1.4653e-02.
x[9]=2.83, y=5.241622. exact=5.223013, error=1.8609e-02.
x[10]=3.14, y=0.021620. exact=-0.000000, error=2.1620e-02.

2、差分法

2.1 狄利克雷边值问题

2.1.1 计算思路

首先考虑形如式（5）的两点狄利克雷边值问题：

$\left\{\begin{matrix} y^{''}(x)+q(x)y(x)=f(x),\space\space\space\space a<x<b,\space\space (5)\\ y(a)=\alpha,\space\space\space\space y(b)=\beta \end{matrix}\right.$

其中 $q(x),f(x)$ 为已知函数且在[a,b]内 $q\leqslant 0$ ， $y(x)$ 为待求函数， $\alpha,\beta$ 为常数。首先对求解域[a,b]等距剖分，得到m+1个节点 $x_{i}=a+ih,i=0,1,\cdot \cdot \cdot ,m$ ，步长为 $h=\frac{b-a}{m}$ 。然后将院复查弱化，使其在内部节点成立，有：

$y^{''}(x_{i})+q(x_{i})y(x_{i})=f(x_{i}),i=1,2,\cdot \cdot \cdot ,m-1$

再对二阶导数值 $y^{''}(x_{i})$ 进行处理，可得到在节点 $x_{i}$ 处的离散方程：

$\frac{y(x_{i+1})-2y(x_{i})+y(x_{i-1})}{h^{2}}-\frac{h^{2}}{12}y^{(4)}(\xi _{i})+q(x_{i})y(x_{i})=f(x_{i}).i=1,2,\cdot \cdot \cdot,m-1$

其中 $\xi_{i}\in (x_{i-1},x_{i+1})$ 。忽略其中的高阶项，用数值解 $y_{i}$ 作为精确解 $y(x_{i})$ 的近似值，结合边界条件，可得：

$\left\{\begin{matrix} \frac{y_{i+1}-2y_{i}+y_{i-1}}{h^{2}}+q(x_{i})y_{i}=f(x_{i}),\space\space i=1,2,\cdot\cdot\cdot,m-1,\\ y_{0}=\alpha,\space\space\space\space y_{m}=\beta \end{matrix}\right.$

这是一个由m-1个未知量 $y_{i}(i=1,2,\cdot\cdot\cdot,m-1)$ 、 $(m-1)$ 个方程构成的线性方程组，可以写成矩阵的形式：

$\begin{pmatrix} -2+h^{2}q(x_{i}) & 1 & & 0 & \\ 1 & -2+h^{2}q(x_{2}) & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2+h^{2}q(x_{m-2}) & 1\\ & 0 & & 1 & -2+h^{2}q(x_{m-1}) \end{pmatrix}$ $\begin{pmatrix} y_{1}\\ y_{2}\\ \vdots \\ y_{m-2}\\ y_{m-1} \end{pmatrix}$ $=$ $\begin{pmatrix} h^{2}f(x_{1})-\alpha\\ h^{2}f(x_{2})\\ \vdots\\ h^{2}f(x_{m-2})\\ h^{2}f(x_{m-1})-\beta \end{pmatrix}$

由于 $q(x)\leqslant 0$ ，所以上述方程组系数矩阵是一个对角占优的三对角矩阵，可使用追赶法直接求解。

2.1.2 算例实现

求解两点边值问题：

$\left\{\begin{matrix} y^{''}(x)-(x-\frac{1}{2})^{2}y(x)=-(x^{2}-x+\frac{5}{4})sinx,\space\space\space\space 0<x<\frac{\pi}{2},\\ y(0)=0,\space\space\space\space y(\frac{\pi}{2})=1 \end{matrix}\right.$

步长h=π/16，计算节点 $x=0,\frac{\pi}{8},\frac{\pi}{4},\frac{3\pi}{8},\frac{\pi}{2}$ 处的数值解，并给出与精确解 $y(x)=sinx$ 的误差。

代码如下：


#include <cmath>
#include <stdlib.h>
#include <stdio.h>
 
int main(int argc, char* argv[])
{
        int i,j,N;
        double a,b,h,Pi,alpha,beta;
        double *matrix_a_array,*matrix_b_array,*x,*y,*rhs,*ans;
        double f(double x);
        double q(double x);
        double *chase_algorithm(double *a, double *b, double *c, double *d, int n);
        double exact(double x);
 
        N=8;
        Pi=3.14159265359;
        a=0.0;
        b=Pi/2.0;
        h=(b-a)/N;
        alpha=0.0;
        beta=1.0;
 
        x=(double*)malloc(sizeof(double)*(N+1));
        for(i=0;i<=N;i++)
                x[i]=a+i*h;
 
        rhs=(double*)malloc(sizeof(double)*(N-1));
        for(i=1;i<N;i++)
                rhs[i-1]=h*h*f(x[i]);
 
        rhs[0]=rhs[0]-alpha;
        rhs[N-2]=rhs[N-2]-beta;
 
        matrix_a_array=(double*)malloc(sizeof(double)*(N-1));
        matrix_b_array=(double*)malloc(sizeof(double)*(N-1));
        for(i=0;i<N-1;i++)
        {
                matrix_a_array[i]=1.0;
                matrix_b_array[i]=h*h*q(x[i+1])-2;
        }
        ans=chase_algorithm(matrix_a_array,matrix_b_array,matrix_a_array,rhs,N-1);
        free(matrix_a_array);free(matrix_b_array);free(rhs);
 
        y=(double*)malloc(sizeof(double)*(N+1));
        y[0]=alpha;
        for(i=1;i<N;i++)
                y[i]=ans[i-1];
        free(ans);
        y[N]=beta;
 
        j=N/4;
        for(i=0;i<=N;i=i+j)
                printf("x=%.2f,y=%f,exact=%f,err=%.4e.\n",x[i],y[i],exact(x[i]),fabs(exact(x[i])-y[i]));
 
        return 0;
}
 
 
double f(double x)
{
        return -(x*x-x+5.0/4.0)*sin(x);
}
double q(double x)
{
        double z;
        z=x-0.5;
        return -z*z;
}
double* chase_algorithm(double *a, double *b, double *c, double *d, int n)
{
        int i;
        double *ans,*g,*w,p;
 
        ans=(double*)malloc(sizeof(double)*n);
        g=(double*)malloc(sizeof(double)*n);
        w=(double*)malloc(sizeof(double)*n);
        g[0]=d[0]/b[0];
        w[0]=c[0]/b[0];
 
        for(i=1;i<n;i++)
        {
                p=b[i]-a[i]*w[i-1];
                g[i]=(d[i]-a[i]*g[i-1])/p;
                w[i]=c[i]/p;
        }
        ans[n-1]=g[n-1];
        i=n-2;
        do
        {
                ans[i]=g[i]-w[i]*ans[i+1];
                i=i-1;
        }while(i>=0);
        free(g);free(w);
 
        return ans;
}
double exact(double x)
{
        return sin(x);
}

N=8时，运行结果如下：


x=0.00,y=0.000000,exact=0.000000,err=0.0000e+00.
x=0.39,y=0.383096,exact=0.382683,err=4.1273e-04.
x=0.79,y=0.707746,exact=0.707107,err=6.3923e-04.
x=1.18,y=0.924409,exact=0.923880,err=5.2906e-04.
x=1.57,y=1.000000,exact=1.000000,err=0.0000e+00.

2.2 混合边值问题

2.2.1 计算思路

考虑形如式（6）的两点混合边值问题：

$\left\{\begin{matrix} y^{''}(x)+q(x)y(x)=f(x),\space\space\space\space a<x<b,\space\space (6)\\ y^{'}(a)+\lambda y(a)=\alpha, \space\space y^{'}(b)+\mu y(b)=\beta \end{matrix}\right.$

式（6）在离散节点处的方程为：

$\left\{\begin{matrix} y^{''}(x_i)+q(x_{i})y(x_{i})=f(x_{i}),\space\space\space\space 1\leqslant i\leqslant m-1,\\ y^{'}(x_{0})+\lambda y(x_{0})=\alpha,\space\space y^{'}(x_{m})+\mu y(x_{m})=\beta \end{matrix}\right.$

对其中内部节点的二阶导数采用二阶差商，起点、终点处的一阶导数分别采用向前、向后差商，忽略高阶项，并用数值解代替精确解，可得差分格式：

$\left\{\begin{matrix} \frac{y_{i+1}-2y_{i}+y_{i-1}}{h}+q(x_{i})y_{i}=f(x_{i}),\space\space i=1,2,\cdot\cdot\cdot,m-1,\\ \frac{y_{1}-y_{0}}{h}+\lambda y_{0}=\alpha,\space\space\space\space \frac{y_{m}-y_{m-1}}{h}+\mu y_{m}=\beta \end{matrix}\right.$

可将上式写成矩阵形式：

$\begin{pmatrix} h\lambda-1& 1 & & 0 & \\ 1 & -2+h^{2}q(x_{1}) & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2+h^{2}q(x_{m-1}) & 1\\ & 0 & & 1 & -(h\mu +1) \end{pmatrix}$ $\begin{pmatrix} y_{0}\\ y_{1}\\ \vdots \\ y_{m-1}\\ y_{m} \end{pmatrix}$ = $\begin{pmatrix} h\alpha\\ h^{2}f(x_{1})\\ \vdots\\ h^{2}f(x_{m-1})\\ -h\beta \end{pmatrix}$

该矩阵为m+1阶的三对角线性方程组，用追赶法求解。

2.2.2 算例实现

求解两点混合边值问题：

$\left\{\begin{matrix} y^{''}(x)-(1+sinx)y(x)=-e^{x}sinx,\space\space\space\space 0<x<1,\\ y^{'}(0)-y(0)=0,\space\space y^{'}(1)+2y(1)=3e \end{matrix}\right.$

步长h=1/4。计算节点 $x=0,\frac{1}{4},\frac{1}{2},\frac{3}{4},1$ 处的数值解，并给出与精确解 $u(x)=e^{x}$ 的误差。

代码如下：


#include <cmath>
#include <stdlib.h>
#include <stdio.h>
 
int main(int argc, char* argv[])
{
        int i,j,N;
        double a,b,h,alpha,beta,lambda,mu;
        double *matrix_a_array,*matrix_b_array,*x,*y,*rhs;
        double f(double x);
        double q(double x);
        double *chase_algorithm(double *a, double *b, double *c, double *d, int n);
        double exact(double x);
 
        a=0.0;
        b=1.0;
        N=4;
        h=(b-a)/N;
        alpha=0.0;
        beta=3*exp(1.0);
        lambda=-1.0;
        mu=2.0;
 
        x=(double*)malloc(sizeof(double)*(N+1));
        for(i=0;i<=N;i++)
                x[i]=a+i*h;
 
        rhs=(double*)malloc(sizeof(double)*(N+1));
        rhs[0]=h*alpha;
        for(i=1;i<N;i++)
                rhs[i]=h*h*f(x[i]);
        rhs[N]=-h*beta;
 
        matrix_a_array=(double*)malloc(sizeof(double)*(N+1));
        matrix_b_array=(double*)malloc(sizeof(double)*(N+1));
        for(i=0;i<=N;i++)
        {
                matrix_a_array[i]=1.0;
                matrix_b_array[i]=h*h*q(x[i])-2;
        }
        matrix_b_array[0]=h*lambda-1.0;
        matrix_b_array[N]=-(h*mu+1.0);
 
        y=chase_algorithm(matrix_a_array,matrix_b_array,matrix_a_array,rhs,N+1);
        free(matrix_a_array);free(matrix_b_array);free(rhs);
 
        j=N/4;
        for(i=0;i<=N;i=i+j)
                printf("x=%.2f,y=%f,exact=%f,err=%f.\n",x[i],y[i],exact(x[i]),fabs(exact(x[i])-y[i]));
        free(x);free(y);
 
        return 0;
}
 
 
double f(double x)
{
        return -exp(x)*sin(x);
}
double q(double x)
{
        return -(1+sin(x));
}
double *chase_algorithm(double *a, double *b, double *c, double *d, int n)
{
        int i;
        double *ans,*g,*w,p;
 
        ans=(double*)malloc(sizeof(double)*n);
        g=(double*)malloc(sizeof(double)*n);
        w=(double*)malloc(sizeof(double)*n);
        g[0]=d[0]/b[0];
        w[0]=c[0]/b[0];
 
        for(i=1;i<n;i++)
        {
                p=b[i]-a[i]*w[i-1];
                g[i]=(d[i]-a[i]*g[i-1])/p;
                w[i]=c[i]/p;
        }
 
        ans[n-1]=g[n-1];
        i=n-2;
        do
        {
                ans[i]=g[i]-w[i]*ans[i+1];
                i=i-1;
        }while(i>=0);
 
        free(g);free(w);
        return ans;
}
double exact(double x)
{
        return exp(x);
}

N=4时，运行结果如下：


x=0.00,y=1.104528,exact=1.000000,err=0.104528.
x=0.25,y=1.380660,exact=1.284025,err=0.096635.
x=0.50,y=1.744578,exact=1.648721,err=0.095856.
x=0.75,y=2.220404,exact=2.117000,err=0.103404.
x=1.00,y=2.839410,exact=2.718282,err=0.121128.

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/Monodyee/article/detail/709592