偏微分方程算法之双曲型方程差分格式编程示例一

猿核试Bug愁

已于 2024-05-20 15:48:51 修改

阅读量596

点赞数 13

分类专栏：偏微分方程（双曲型）文章标签：算法

于 2024-05-11 11:54:02 首次发布

本文链接：https://blog.csdn.net/L_peanut/article/details/138668938

版权

偏微分方程（双曲型）专栏收录该内容

11 篇文章 12 订阅 ¥9.90 ¥99.00

订阅专栏

超级会员免费看

该博客详细介绍了双曲型偏微分方程的显格式、隐格式和紧差分格式，通过C++编程示例展示了如何对这三种格式进行求解，并给出了不同步长下的计算结果，帮助读者深入理解差分格式。

摘要由CSDN通过智能技术生成

本专栏对双曲型偏微分方程的显格式、隐格式、紧差分格式方法进行了介绍，并给出相应格式的理论推导过程。为加深对差分格式的理解，分别对三种方法进行C++编程示例。

采用不同的差分格式对同一二阶双曲型方程初边值问题进行编程求解。问题如下：

$\left\{\begin{matrix} \frac{\partial^{2}u}{\partial t^{2}}-\frac{\partial^{2}u}{\partial x^{2}}=\frac{2t(1-2x^{2}-3x^{4})}{(1+x^{2})^{4}},-1<x<1,0<t\leqslant 2,\\ u(x,0)=0,\frac{\partial u}{\partial t}(x,0)=\frac{1}{1+x^{2}},-1\leqslant x\leqslant 1,\\ u(-1,t)=u(1,t)=\frac{t}{2},0<t\leqslant 2 \end{matrix}\right.$

已知问题的精确解为 $u(x,t)=\frac{t}{1+x^{2}}$ 。

一、显格式

分别取步长 $\tau_{1}=\frac{1}{50},h_{1}=\frac{1}{20}$ 和 $\tau_{1}=\frac{1}{100},h_{1}=\frac{1}{40}$ ，给出在节点 $(0.2,1+0.2i),i=0,1,\cdot\cdot\cdot,5$ 处的数值解及误差。

1.1 C++代码

#include <cmath>
#include <stdlib.h>
#include <stdio.h>

int main(int argc, char *argv[])
{
        int m,n,i,j,k;
        double a,r,h,tau,*x,*t,**u;
        double phi(double x);
        double psi(double x);
        double alpha(double t);
        double beta(double t);
        double f(double x, double t);
        double exact(double x, double t);

        m=20;
        n=50;
        a=1.0;
        h=2.0/m;
        tau=2.0/n;
        r=a*tau/h;
        r=r*r;
        printf("r=%.4f.\n",r);
        printf("h=%f, tau=%f.\n",h,tau);

        x=(double*)malloc(sizeof(double)*(m+1));
        for(i=0;i<=m;i++)
                x[i]=-1.0+i*h;

        t=(double*)malloc(sizeof(double)*(n+1));
        for(k=0;k<=n;k++)
                t[k]=k*tau;

        u=(double**)malloc(sizeof(double*)*(m+1));
        for(i=0;i<=m;i++)
                u[i]=(double*)malloc(sizeof(double)*(n+1));

        for(i=0;i<=m;i++)
                u[i][0]=phi(x[i]);

        for(k=1;k<=n;k++)
        {
                u[0][k]=alpha(t[k]);
                u[m][k]=beta(t[k]);
        }

        for(i=1;i<m;i++)
                u[i][1]=(r*u[i-1][0]+2*(1-r)*u[i][0]+r*u[i+1][0]+tau*tau*f(x[i],t[0])+2*tau*psi(x[i]))/2.0;

        for(k=1;k<n;k++)
        {
                for(i=1;i<m;i++)
                {
                        u[i][k+1]=r*u[i-1][k]+2*(1-r)*u[i][k]+r*u[i+1][k]-u[i][k-1]+tau*tau*f(x[i],t[k]);
                }
        }

        j=3*m/5;
        k=n/10;
        for(i=n/2;i<=n;i=i+k)
                printf("(x,t)=(%.2f,%.2f), y=%f, error=%.4e.\n",x[j],t[i],u[j][i],fabs(u[j][i]-exact(x[j],t[i])));

        for(i=0;i<=m;i++)
                free(u[i]);
        free(u);free(x);free(t);

        return 0;
}




double phi(double x)
{
        return 0.0;
}
double psi(double x)
{
        return 1.0/(1.0+x*x);
}
double alpha(double t)
{
        return t/2.0;
}
double beta(double t)
{
        return t/2.0;
}
double f(double x, double t)
{
        return 2*t*(1.0-2*x*x-3.0*pow(x,4))/pow((1.0+x*x),4);
}
double exact(double x, double t)
{
        return t/(1.0+x*x);
}

1.2 计算结果

当 $\tau_{1}=\frac{1}{50},h_{1}=\frac{1}{20}$ 时，计算结果如下：

r=0.1600.
h=0.100000, tau=0.040000.
(x,t)=(0.20,1.00), y=0.962434, error=8.9599e-04.
(x,t)=(0.20,1.20), y=1.155106, error=1.2594e-03.
(x,t)=(0.20,1.40), y=1.347825, error=1.6715e-03.
(x,t)=(0.20,1.60), y=1.540629, error=2.1678e-03.
(x,t)=(0.20,1.80), y=1.733539, error=2.7700e-03.
(x,t)=(0.20,2.00), y=1.926523, error=3.4464e-03.

当 $\tau_{1}=\frac{1}{100},h_{1}=\frac{1}{40}$ 时，计算结果如下：

r=0.1600.
h=0.050000, tau=0.020000.
(x,t)=(0.20,1.00), y=0.961762, error=2.2313e-04.
(x,t)=(0.20,1.20), y=1.154160, error=3.1358e-04.
(x,t)=(0.20,1.40), y=1.346570, error=4.1629e-04.
(x,t)=(0.20,1.60), y=1.539002, error=5.4082e-04.
(x,t)=(0.20,1.80), y=1.731461, error=6.9173e-04.
(x,t)=(0.20,2.00), y=1.923938, error=8.6061e-04.

二、隐格式

分别取步长 $\tau_{1}=\frac{1}{50},h_{1}=\frac{1}{50}$ 和 $\tau_{1}=\frac{1}{100},h_{1}=\frac{1}{100}$ ，给出在节点 $(0.2,1+0.2i),i=0,1,\cdot\cdot\cdot,5$ 处的数值解及误差。

2.1 C++代码

#include <cmath>
#include <stdlib.h>
#include <stdio.h>

int main(int argc, char *argv[])
{
        int m,n,i,j,k;
        double a,h,r,tau,*x,*t,**u,*a1,*b,*c,*d,*ans;
        double phi(double x);
        double psi(double x);
        double alpha(double t);
        double beta(double t);
        double f(double x, double t);
        double exact(double x, double t);
        double *chase_algorithm(double *a, double *b, double *c, double *d, int n);

        m=50;
        n=50;
        a=1.0;
        h=2.0/m;
        tau=2.0/n;
        r=a*tau/h;
        r=r*r;
        printf("r=%.4f.\n",r);
        printf("h=%f,tau=%f.\n",h,tau);

        x=(double*)malloc(sizeof(double)*(m+1));
        for(i=0;i<=m;i++)
                x[i]=-1.0+i*h;

        t=(double*)malloc(sizeof(double)*(n+1));
        for(k=0;k<=n;k++)
                t[k]=k*tau;

        u=(double**)malloc(sizeof(double*)*(m+1));
        for(i=0;i<=m;i++)
                u[i]=(double*)malloc(sizeof(double)*(n+1));

        for(i=0;i<=m;i++)
                u[i][0]=phi(x[i]);

        for(k=1;k<=n;k++)
        {
                u[0][k]=alpha(t[k]);
                u[m][k]=beta(t[k]);
        }

        for(i=1;i<m;i++)
                u[i][1]=(r*u[i-1][0]+2*(1-r)*u[i][0]+r*u[i+1][0]+tau*tau*f(x[i],t[0])+2*tau*psi(x[i]))/2.0;

        a1=(double *)malloc(sizeof(double)*(m-1));
        b=(double *)malloc(sizeof(double)*(m-1));
        c=(double *)malloc(sizeof(double)*(m-1));
        d=(double *)malloc(sizeof(double)*(m-1));
        ans=(double *)malloc(sizeof(double)*(m-1));

        for(k=1;k<n;k++)
        {
                for(i=1;i<m;i++)
                {
                        d[i-1]=r*(u[i-1][k-1]+u[i+1][k-1])/2.0-(1+r)*u[i][k-1]+2*u[i][k]+tau*tau*f(x[i],t[k]);
                        a1[i-1]=-0.5*r;
                        b[i-1]=1.0+r;
                        c[i-1]=a1[i-1];
                }
                d[0]=d[0]+0.5*r*u[0][k+1];
                d[m-2]=d[m-2]+0.5*r*u[m][k+1];
                ans=chase_algorithm(a1,b,c,d,m-1);
                for(i=0;i<m-1;i++)
                        u[i+1][k+1]=ans[i];
        }
        free(ans);

        j=3*m/5;
        k=n/10;
        for(i=n/2;i<=n;i=i+k)
                printf("(x,t)=(%.2f,%.2f), y=%f, error=%.4e.\n",x[j],t[i],u[j][i],fabs(u[j][i]-exact(x[j],t[i])));

        for(i=0;i<=m;i++)
                free(u[i]);
        free(u);free(a1);free(b);free(c);free(d);free(x);free(t);

        return 0;
}



double phi(double x)
{
        return 0.0;
}
double psi(double x)
{
        return 1.0/(1.0+x*x);
}
double alpha(double t)
{
        return t/2.0;
}
double beta(double t)
{
        return t/2.0;
}
double f(double x, double t)
{
        return 2*t*(1.0-2*x*x-3.0*pow(x,4))/pow((1.0+x*x),4);
}
double exact(double x, double t)
{
        return t/(1.0+x*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;
}

2.2 计算结果

当 $\tau_{1}=\frac{1}{50},h_{1}=\frac{1}{50}$ 时，计算结果如下：

r=1.0000.
h=0.040000,tau=0.040000.
(x,t)=(0.20,1.00), y=0.961681, error=1.4270e-04.
(x,t)=(0.20,1.20), y=1.154047, error=2.0072e-04.
(x,t)=(0.20,1.40), y=1.346420, error=2.6654e-04.
(x,t)=(0.20,1.60), y=1.538807, error=3.4573e-04.
(x,t)=(0.20,1.80), y=1.731211, error=4.4173e-04.
(x,t)=(0.20,2.00), y=1.923627, error=5.4959e-04.

当 $\tau_{1}=\frac{1}{100},h_{1}=\frac{1}{100}$ 时，计算结果如下：

r=1.0000.
h=0.020000,tau=0.020000.
(x,t)=(0.20,1.00), y=0.961574, error=3.5660e-05.
(x,t)=(0.20,1.20), y=1.153896, error=5.0125e-05.
(x,t)=(0.20,1.40), y=1.346220, error=6.6550e-05.
(x,t)=(0.20,1.60), y=1.538548, error=8.6461e-05.
(x,t)=(0.20,1.80), y=1.730880, error=1.1058e-04.
(x,t)=(0.20,2.00), y=1.923215, error=1.3758e-04.

三、紧差分格式

分别取步长 $\tau_{1}=\frac{1}{20},h_{1}=\frac{1}{20}$ 和 $\tau_{1}=\frac{1}{80},h_{1}=\frac{1}{40}$ ，给出在节点 $(0.2,1+0.2i),i=0,1,\cdot\cdot\cdot,5$ 处的数值解及误差。

3.1 C++代码

#include <cmath>
#include <stdlib.h>
#include <stdio.h>

int main(int argc, char* argv[])
{
        int m,n,i,j,k;
        double a,h,r,tau,c1,c2,*x,*t,**u,*a1,*b,*c,*d,*ans;
        double phi(double x);
        double psi(double x);
        double alpha(double t);
        double beta(double t);
        double f(double x, double t);
        double exact(double x, double t);
        double *chase_algorithm(double *a, double *b, double *c, double *d, int n);

        m=20;
        n=20;
        a=1.0;
        h=2.0/m;
        tau=2.0/n;
        r=a*tau/h;
        r=r*r;
        printf("r=%.4f.\n",r);
        printf("h=%f,tau=%f.\n",h,tau);

        x=(double*)malloc(sizeof(double)*(m+1));
        for(i=0;i<=m;i++)
                x[i]=-1.0+i*h;

        t=(double*)malloc(sizeof(double)*(n+1));
        for(k=0;k<=n;k++)
                t[k]=k*tau;

        u=(double**)malloc(sizeof(double*)*(m+1));
        for(i=0;i<=m;i++)
                u[i]=(double*)malloc(sizeof(double)*(n+1));

        for(i=0;i<=m;i++)
                u[i][0]=phi(x[i]);

        for(k=1;k<=n;k++)
        {
                u[0][k]=alpha(t[k]);
                u[m][k]=beta(t[k]);
        }

        for(i=1;i<m;i++)
                u[i][1]=(r*u[i-1][0]+2*(1-r)*u[i][0]+r*u[i+1][0]+tau*tau*f(x[i],t[0])+2*tau*psi(x[i]))/2.0;

        a1=(double *)malloc(sizeof(double)*(m-1));
        b=(double *)malloc(sizeof(double)*(m-1));
        c=(double *)malloc(sizeof(double)*(m-1));
        d=(double *)malloc(sizeof(double)*(m-1));
        ans=(double *)malloc(sizeof(double)*(m-1));
        c1=1.0-6*r;
        c2=10.0+12*r;

        for(k=1;k<n;k++)
        {
                for(i=1;i<m;i++)
                {
                        d[i-1]=(-c1)*(u[i-1][k-1]+u[i+1][k-1])-c2*u[i][k-1]+2*(u[i-1][k]+10*u[i][k]+u[i+1][k])+tau*tau*(f(x[i-1],t[k])+10*f(x[i],t[k])+f(x[i+1],t[k]));
                        a1[i-1]=c1;
                        b[i-1]=c2;
                        c[i-1]=a1[i-1];
                }
                d[0]=d[0]-c1*u[0][k+1];
                d[m-2]=d[m-2]-c1*u[m][k+1];
                ans=chase_algorithm(a1,b,c,d,m-1);
                for(i=0;i<m-1;i++)
                        u[i+1][k+1]=ans[i];
        }
        free(ans);

        j=3*m/5;
        k=n/10;
        for(i=n/2;i<=n;i=i+k)
                printf("(x,t)=(%.2f,%.2f), y=%f, error=%.4e.\n",x[j],t[i],u[j][i],fabs(u[j][i]-exact(x[j],t[i])));

        for(i=0;i<=m;i++)
                free(u[i]);
        free(u);free(a1);free(b);free(c);free(d);free(x);free(t);

        return 0;
}


double phi(double x)
{
        return 0.0;
}
double psi(double x)
{
        return 1.0/(1.0+x*x);
}
double alpha(double t)
{
        return t/2.0;
}
double beta(double t)
{
        return t/2.0;
}
double f(double x, double t)
{
        return 2*t*(1.0-2*x*x-3.0*pow(x,4))/pow((1.0+x*x),4);
}
double exact(double x, double t)
{
        return t/(1.0+x*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;
}

3.1 计算结果

当 $\tau_{1}=\frac{1}{20},h_{1}=\frac{1}{20}$ 时，计算结果如下：

r=1.0000.
h=0.100000,tau=0.100000.
(x,t)=(0.20,1.00), y=0.961543, error=4.9097e-06.
(x,t)=(0.20,1.20), y=1.153852, error=6.0706e-06.
(x,t)=(0.20,1.40), y=1.346161, error=7.0463e-06.
(x,t)=(0.20,1.60), y=1.538470, error=8.1881e-06.
(x,t)=(0.20,1.80), y=1.730779, error=9.8633e-06.
(x,t)=(0.20,2.00), y=1.923089, error=1.2176e-05.

当 $\tau_{1}=\frac{1}{80},h_{1}=\frac{1}{40}$ 时，计算结果如下:

r=0.2500.
h=0.050000,tau=0.025000.
(x,t)=(0.20,1.00), y=0.961539, error=3.0354e-07.
(x,t)=(0.20,1.20), y=1.153847, error=3.7180e-07.
(x,t)=(0.20,1.40), y=1.346154, error=4.3037e-07.
(x,t)=(0.20,1.60), y=1.538462, error=5.0731e-07.
(x,t)=(0.20,1.80), y=1.730770, error=6.2824e-07.
(x,t)=(0.20,2.00), y=1.923078, error=7.8048e-07.