简单的差分格式解一个一元二阶常微分方程_差分求解二阶微分方程-CSDN博客

本文链接：https://blog.csdn.net/2301_76273740/article/details/140226328

对于一个如下的一元二阶常微分方程
$\begin{cases} -u''+q\left( x \right) u=f\left( x \right) \,\, a<x<b\\ u\left( a \right) =\alpha \,\,u\left( b \right) =\beta\\ \end{cases}$
其中二阶导数可以换成差分格式
$u''\left( x_i \right) =\frac{1}{h^2}\left( u\left( x_{i-1} \right) -2u\left( x_i \right) +u\left( x_{i+1} \right) \right)$
将 $u(x_i)$ 简写为 $u_i$
$-\frac{u_{i-1}}{h^2}+\left( q\left( x_i \right) +\frac{2}{h^2} \right) u_i-\frac{u_{i+1}}{h^2}=f\left( x_i \right)$
将a到b的区域平均分成m份，解以上的式子可以看做解一个三对角方程组
$\left( \begin{matrix} \begin{matrix} 2+h^2q(x_1)& -1& \\ -1& 2+h^2q(x_2)& -1\\ & \ddots& \ddots\\ \end{matrix}& \begin{matrix} & & \\ & & \\ \ddots& & \\ \end{matrix}\\ \begin{matrix} & & \\ & & \\ & & \\ \end{matrix}& \begin{matrix} -1& 2+h^2q(x_{m-2})& -1\\ & -1& 2+h^2q(x_{m-1})\\ & & \\ \end{matrix}\\ \end{matrix} \right) \left( \begin{array}{c} u_1\\ u_2\\ \begin{array}{c} \vdots\\ u_{m-2}\\ \end{array}\\ u_{m-1}\\ \end{array} \right) =\left( \begin{array}{c} h^2f\left( x_1 \right) +\alpha\\ h^2f\left( x_2 \right)\\ \begin{array}{c} \vdots\\ h^2f\left( x_{m-2} \right)\\ \end{array}\\ h^2f\left( x_{m-1} \right) +\beta\\ \end{array} \right)$
这里是左右两边同时乘以 $h^2$ 之后展开成的矩阵形式，取 $f(x)=e^x(\sin(x)-2\cos(x))$ $q (x) = 1$ 代码如下


#include<stdio.h>
#include<math.h>
#define M_PI 3.14159265358979323846

#define M 20 //m等分
#define N M-1
/*
a、b、c是三对角矩阵的下对角线、主对角线和上对角线的元素，d是右侧向量，x是结果向量
*/
void thomas(double a[], double b[], double c[], double d[], double x[])
{
    double y[N],p[N],q[N];
    int i;
    p[0]=b[0];
    for(i=0;i<N-1;i++){
        q[i]=c[i]/p[i];
        p[i+1]=b[i+1]-a[i]*q[i];
    }
    y[0]=d[0]/p[0];
    for(i=1;i<N;i++){
        y[i]=(d[i]-a[i-1]*y[i-1])/p[i];
    }
    x[N-1]=y[N-1];
    for(i=N-2;i>=0;i--){
        x[i]=y[i]-q[i]*x[i+1];
    }
}

double f(double x){
    return exp(x)*(sin(x)-2*cos(x));
}

double q(double x){
    return 1;
}
int main(){
    double a=0.0;//下限
    double b=M_PI;//上限
    double h=(b-a)/M;//步长
    double alpha=0.0;
    double beta=0.0;
    double zhu[N],xia[N-1],shang[N-1],you[N],x[N];
    int i;
    for(i=0;i<N;i++){
        zhu[i]=q(a+i*h+h)*h*h+2;
        you[i]=f(a+i*h+h)*h*h;
    }
    you[0]+=alpha;
    you[N-1]+=beta;
    for(i=0;i<N-1;i++){
        xia[i]=-1;
        shang[i]=-1;
    }
    thomas(xia,zhu,shang,you,x);
    for(i=0;i<N;i++){
        printf("u[%d]=%lf\n",i+1,x[i]);
    }
    return 0;
}