偏微分方程算法之二阶双曲型方程紧差分方法

2024-04-22 16:36

本文主要是介绍偏微分方程算法之二阶双曲型方程紧差分方法,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

目录

一、研究目标

二、理论推导

三、算例实现


一、研究目标

        前面我们已经介绍了二阶双曲型方程显式、隐式差分格式,可否像抛物型方程一样,构建更高精度的差分格式。接下来我们介绍紧差分格式。这里继续以非齐次二阶双曲型偏微分方程的初边值问题为研究对象:

\left\{\begin{matrix} \frac{\partial^{2}u(x,t)}{\partial t^{2}}-a^{2}\frac{\partial^{2}u(x,t)}{\partial x^{2}}=f(x,t),0<x<1,0<t\leqslant T,\\ u(x,0)=\varphi(x),\frac{\partial u}{\partial t}(x,0)=\Psi(x),0\leqslant x\leqslant 1,\space\space(1)\\ u(0,t)=\alpha(t),u(1,t)=\beta(t),0<t\leqslant T \end{matrix}\right.

公式(1)中u表示一个与时间t和位置x有关的待求波函数,\varphi(x),\Psi(x),\alpha(t),\beta(t)及方程右端项函数f(x,t)都是已知函数,a,T是非零常数。

二、理论推导

        第一步:网格剖分。对矩形求解域0\leqslant x\leqslant 1,0\leqslant t\leqslant T进行等距剖分,即

x_{j}=jh(j=0,1,\cdot\cdot\cdot,m),t_{k}=k\tau(k=0,1,\cdot\cdot\cdot,n)

        第二步:弱化原方程。将原来的连续方程离散到网格节点上成立,得到离散方程:

\left\{\begin{matrix} \frac{\partial^{2}u}{\partial t^{2}}|_{(x_{j},t_{k})}-a^{2}\frac{\partial ^{2}u}{\partial x^{2}}|_{(x_{j},t_{k})}=f(x_{j},t_{k}),0<j<m,0<k\leqslant n,\\ u(x_{j},t_{0})=\varphi(x_{j}),\frac{\partial u}{\partial t}(x_{j},t_{0})=\Psi(x_{j}),0\leqslant j\leqslant m,\space\space(2)\\ u(x_{0},t_{k})=\alpha(t_{k}),u(x_{m},t_{k})=\beta(t_{k}),1\leqslant k\leqslant n \end{matrix}\right.

        第三步:对偏导数进行更高精度近似。由泰勒公式(固定时间t不变):

u(x_{j-1},t_{k})=(u-h\frac{\partial u}{\partial x}+\frac{h^{2}}{2}\frac{\partial ^{2}u}{\partial x^{2}}-\frac{h^{3}}{6}\frac{\partial^{3}u}{\partial x^{3}}+\frac{h^{4}}{24}\frac{\partial^{4}u}{\partial x^{4}}-\frac{h^{5}}{120}\frac{\partial^{5}u}{\partial x^{5}})|_{(x_{j},t_{k})}+O(h^{6})

u(x_{j+1},t_{k})=(u+h\frac{\partial u}{\partial x}+\frac{h^{2}}{2}\frac{\partial ^{2}u}{\partial x^{2}}+\frac{h^{3}}{6}\frac{\partial^{3}u}{\partial x^{3}}+\frac{h^{4}}{24}\frac{\partial^{4}u}{\partial x^{4}}+\frac{h^{5}}{120}\frac{\partial^{5}u}{\partial x^{5}})|_{(x_{j},t_{k})}+O(h^{6})

将上面两式相加可得

u(x_{j+1},t_{k})-2u(x_{j},t_{k})+u(x_{j-1},t_{k})=h^{2}\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k})}+\frac{h^{4}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k})}+O(h^{6})

有           \frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k})}=\frac{u(x_{j+1},t_{k})-2u(x_{j},t_{k})+u(x_{j-1},t_{k})}{h^{2}}-\frac{h^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k})}+O(h^{4})

类似的有

\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k-1})}=\frac{u(x_{j+1},t_{k-1})-2u(x_{j},t_{k-1})+u(x_{j-1},t_{k-1})}{h^{2}}-\frac{h^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k-1})}+O(h^{4})

\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k+1})}=\frac{u(x_{j+1},t_{k+1})-2u(x_{j},t_{k+1})+u(x_{j-1},t_{k+1})}{h^{2}}-\frac{h^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k+1})}+O(h^{4})

将上面两式相加后除以2得

\frac{1}{2}(\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k-1})}+\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k+1})})=\frac{u(x_{j+1},t_{k-1})-2u(x_{j},t_{k-1})+u(x_{j-1},t_{k+1})}{2h^{2}}+\frac{u(x_{j+1},t_{k+1})-2u(x_{j},t_{k+1})+u(x_{j-1},t_{k+1})}{2h{^{2}}}-\frac{h^{2}}{12}\frac{1}{2}(\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k-1})}+\frac{\partial ^{4}u}{\partial x^{4}}|_{(x_{j},t_{k+1})})+O(h^{4})

从而         \frac{\partial^{2}u}{\partial x^{2}}|_{(x_{j},t_{k})}+O(\tau^{2})=\frac{u(x_{j+1},t_{k-1})-2u(x_{j},t_{k-1})+u(x_{j-1},t_{k-1})}{2h^{2}}+\frac{u(x_{j+1},t_{k+1})-2u(x_{j},t_{_{k+1}})+u(x_{j-1},t_{k+1})}{2h^{2}}-\frac{h^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{j},t_{k})}+O(h^{4}+\tau^{2}h^{2}) \space\space(3)

现令

D=u(x_{j+1},t_{k-1})-2u(x_{j},t_{k-1})+u(x_{j-1},t_{k-1})+u(x_{j+1},t_{k+1})-2u(x_{j},t_{k+1})+u(x_{j-1},t_{k+1})\frac{\partial^{2}u}{\partial x^{2}}=v(x,t)

则公式(3)可写作

v(x_{j},t_{k})=\frac{D}{2h^{2}}-\frac{h^{2}}{12}\frac{\partial^{2}v}{\partial x^{2}}|_{(x_{j},t_{k})}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})

=\frac{D}{2h^{2}}-\frac{h^{2}}{12}\frac{v(x_{j+1},t_{k})-2v(x_{j},t_{k})+v(x_{j-1},t_{k})}{h^{2}}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})

整理可得

\frac{v(x_{j+1},t_{k})+10v(x_{j},t_{k})+v(x_{j-1},t_{k})}{12}=\frac{D}{2h^{2}}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})\space\space\space\space(4)

        由于原双曲型方程为\frac{\partial^{2}u}{\partial t^{2}}-a^{2}v(x,t)=f(x,t),也即v(x,t)=\frac{1}{a^{2}}(\frac{\partial^{2}u}{\partial t^{2}}-f(x,t))

        公式(4)可改写为

\frac{1}{12a^{2}}(\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{j-1},t_{k})}-f(x_{j-1},t_{k})+10\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{j},t_{k})}-10f(x_{j},t_{k})+\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{j+1},t_{k})}-f(x_{j+1},t_{k}))=\frac{D}{2h^{2}}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})

再利用中心差商近似

\frac{1}{12a^{2}}(\frac{u(x_{j-1},t_{k+1})-2u(x_{j-1},t_{k})+u(x_{j-1},t_{k-1})}{\tau^{2}}-f(x_{j-1},t_{k})+10\frac{u(x_{j},t_{k+1})-2u(x_{j},t_{k})+u(x_{j},t_{k-1})}{\tau^{2}}-10f(x_{j},t_{k})+\frac{u(x_{j+1},t_{k+1})-2u(x_{j+1},t_{k}+u(x_{j+1},t_{k-1}))}{\tau^{2}}-f(x_{j+1},t_{k}))=\frac{D}{2h^{2}}+O(\tau^{2}h^{2}+h^{4}+\tau^{2})

        上式中用数值解代替精确解并忽略高阶项,可得

\frac{1}{12a^{2}\tau^{2}}(u^{k+1}_{j-1}-2u^{k}_{j-1}+u^{k+1}_{j-1}+10(u^{k+1}_{j}-2u^{k}_{j}+u^{k-1}_{j})+u^{k+1}_{j+1}-2u^{k}_{j+1}+u^{k-1}_{j+1})=\frac{1}{2}h^{2}(u^{k-1}_{j+1}-2u^{k-1}_{j}+u^{k-1}_{j-1}+u^{k+1}_{j-1}-2u^{k+1}_{j}+u^{k+1}_{j-1})+\frac{1}{12a^{2}}(f^{k}_{j+1}+10f^{k}_{j}+f^{k}_{j-1})

其中,f^{k}_{l}=f(x_{l},t_{k}),l=j-1,j,j+1

        联合初边值条件,可得到以下紧差分格式:

\left\{\begin{matrix} (6r-1)u^{k+1}_{j-1}+(10+12r)u^{k+1}_{j}+(1-6r)u^{k+1}_{j+1}=(6r-1)u^{k-1}_{j-1}-(12r+10)u^{k-1}_{j}+(6r-1)u^{k-1}_{j+1}+\\2(u^{k}_{j-1}+10u^{k}_{j}+u^{k}_{j+1})+\tau^{2}(f^{k}_{j-1}+10f^{k}_{j}+f^{k}_{j+1}),1\leqslant i\leqslant m-1,1\leqslant k\leqslant n-1,\\ u^{0}_{j}=\varphi(_{j}),0\leqslant j\leqslant m,\\ u^{1}_{j}=(ru^{0}_{j-1}+2(1-r)u^{0}_{j}+ru^{0}_{j+1}+\tau^{2}f(x_{j},t_{0})+2\tau\Psi(x_{j}))/2,1\leqslant j\leqslant m-1,\\ u^{k}_{0}=\alpha(t_{k}),u^{k}_{m}=\beta(t_{k}),1\leqslant k\leqslant n \end{matrix}\right.

        其中r=a^{2}\tau^{2}/h^{2}>0,局部截断误差为O(\tau^{2}+h^{4}),关于时间t是二阶的,关于空间x是四阶的。

三、算例实现

        紧差分格式计算双曲型偏微分方程初边值问题:

\left\{\begin{matrix} \frac{\partial^{2}u(x,t)}{\partial t^{2}}-\frac{\partial^{2}u(x,t)}{\partial x^{2}}=2e^{t}sinx,0<x<\pi,0<t\leqslant 1,\\ u(x,0)=sinx,\frac{\partial u}{\partial t}(x,0)=sinx,0\leqslant x\leqslant \pi,\\ u(0,t)=0,u(\pi,t)=0,0<t\leqslant 1 \end{matrix}\right.

已知其精确解为u(x,t)=e^{t}sinx。分别取步长为\tau_{1}=1/50,h_{1}=\pi/200\tau_{1}=1/100,h_{1}=\pi/400,给出节点(\frac{i\pi}{10},\frac{4}{5}),i=1,\cdot\cdot\cdot,5处的数值解和误差。

代码如下:


#include <cmath>
#include <stdlib.h>
#include <stdio.h>int main(int argc, char* argv[])
{int i,j,k,m,n;double a,h,r,tau,pi,c1,c2;double *x,*t,**u,*a1,*b,*c,*d,*ans;double phi(double x);double ddphi(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=200;n=50;a=1.0;pi=3.14159265359;h=pi/m;tau=1.0/n;r=a*tau/h;r=r*r;printf("r=%.4f.\n",r);x=(double*)malloc(sizeof(double)*(m+1));for(i=0;i<=m;i++)x[i]=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);k=4*n/5;j=m/10;for(i=j;i<=m/2;i=i+j)printf("(x,t)=(%.2f,%.2f),y=%f,error=%.4e.\n",x[i],t[k],u[i][k],fabs(u[i][k]-exact(x[i],t[k])));free(a1);free(b);free(c);free(d);free(x);free(t);return 0;
}double phi(double x)
{return sin(x);
}
double psi(double x)
{return sin(x);
}
double alpha(double t)
{return 0.0;
}
double beta(double t)
{return 0.0;
}
double f(double x, double t)
{return 2*sin(x)*exp(t);
}
double exact(double x, double t)
{return sin(x)*exp(t);
}
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;
}

\tau_{1}=1/50,h_{1}=\pi/200时,计算结果如下:

r=1.6211.
(x,t)=(0.31,0.80),y=0.687686,error=4.3538e-05.
(x,t)=(0.63,0.80),y=1.308057,error=8.2815e-05.
(x,t)=(0.94,0.80),y=1.800386,error=1.1398e-04.
(x,t)=(1.26,0.80),y=2.116481,error=1.3400e-04.
(x,t)=(1.57,0.80),y=2.225400,error=1.4089e-04.

\tau_{1}=1/100,h_{1}=\pi/400时,计算结果如下:

r=0.4053.
(x,t)=(0.31,0.80),y=0.687727,error=2.7423e-06.
(x,t)=(0.63,0.80),y=1.308135,error=5.2162e-06.
(x,t)=(0.94,0.80),y=1.800493,error=7.1794e-06.
(x,t)=(1.26,0.80),y=2.116607,error=8.4399e-06.
(x,t)=(1.57,0.80),y=2.225532,error=8.8743e-06.

        从计算结果可知,当时间步长减小为原来的1/4、空间步长减小为原来的1/2时,误差减小为原来的1/16。

这篇关于偏微分方程算法之二阶双曲型方程紧差分方法的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/926292

相关文章

Linux换行符的使用方法详解

《Linux换行符的使用方法详解》本文介绍了Linux中常用的换行符LF及其在文件中的表示,展示了如何使用sed命令替换换行符,并列举了与换行符处理相关的Linux命令,通过代码讲解的非常详细,需要的... 目录简介检测文件中的换行符使用 cat -A 查看换行符使用 od -c 检查字符换行符格式转换将

SpringBoot实现数据库读写分离的3种方法小结

《SpringBoot实现数据库读写分离的3种方法小结》为了提高系统的读写性能和可用性,读写分离是一种经典的数据库架构模式,在SpringBoot应用中,有多种方式可以实现数据库读写分离,本文将介绍三... 目录一、数据库读写分离概述二、方案一:基于AbstractRoutingDataSource实现动态

Java中的String.valueOf()和toString()方法区别小结

《Java中的String.valueOf()和toString()方法区别小结》字符串操作是开发者日常编程任务中不可或缺的一部分,转换为字符串是一种常见需求,其中最常见的就是String.value... 目录String.valueOf()方法方法定义方法实现使用示例使用场景toString()方法方法

Java中List的contains()方法的使用小结

《Java中List的contains()方法的使用小结》List的contains()方法用于检查列表中是否包含指定的元素,借助equals()方法进行判断,下面就来介绍Java中List的c... 目录详细展开1. 方法签名2. 工作原理3. 使用示例4. 注意事项总结结论:List 的 contain

macOS无效Launchpad图标轻松删除的4 种实用方法

《macOS无效Launchpad图标轻松删除的4种实用方法》mac中不在appstore上下载的应用经常在删除后它的图标还残留在launchpad中,并且长按图标也不会出现删除符号,下面解决这个问... 在 MACOS 上,Launchpad(也就是「启动台」)是一个便捷的 App 启动工具。但有时候,应

SpringBoot日志配置SLF4J和Logback的方法实现

《SpringBoot日志配置SLF4J和Logback的方法实现》日志记录是不可或缺的一部分,本文主要介绍了SpringBoot日志配置SLF4J和Logback的方法实现,文中通过示例代码介绍的非... 目录一、前言二、案例一:初识日志三、案例二:使用Lombok输出日志四、案例三:配置Logback一

Python实现无痛修改第三方库源码的方法详解

《Python实现无痛修改第三方库源码的方法详解》很多时候,我们下载的第三方库是不会有需求不满足的情况,但也有极少的情况,第三方库没有兼顾到需求,本文将介绍几个修改源码的操作,大家可以根据需求进行选择... 目录需求不符合模拟示例 1. 修改源文件2. 继承修改3. 猴子补丁4. 追踪局部变量需求不符合很

SpringBoot实现MD5加盐算法的示例代码

《SpringBoot实现MD5加盐算法的示例代码》加盐算法是一种用于增强密码安全性的技术,本文主要介绍了SpringBoot实现MD5加盐算法的示例代码,文中通过示例代码介绍的非常详细,对大家的学习... 目录一、什么是加盐算法二、如何实现加盐算法2.1 加盐算法代码实现2.2 注册页面中进行密码加盐2.

mysql出现ERROR 2003 (HY000): Can‘t connect to MySQL server on ‘localhost‘ (10061)的解决方法

《mysql出现ERROR2003(HY000):Can‘tconnecttoMySQLserveron‘localhost‘(10061)的解决方法》本文主要介绍了mysql出现... 目录前言:第一步:第二步:第三步:总结:前言:当你想通过命令窗口想打开mysql时候发现提http://www.cpp

Mysql删除几亿条数据表中的部分数据的方法实现

《Mysql删除几亿条数据表中的部分数据的方法实现》在MySQL中删除一个大表中的数据时,需要特别注意操作的性能和对系统的影响,本文主要介绍了Mysql删除几亿条数据表中的部分数据的方法实现,具有一定... 目录1、需求2、方案1. 使用 DELETE 语句分批删除2. 使用 INPLACE ALTER T