双曲方程基于matlab的数值解法

双曲型方程基于MATLAB的数值解法（数学1201,陈晓云,41262022）一：一阶双曲型微分方程的初边值问题0,01,01.(,0)cos(),01.(0,)(1,)cos(),01.uuxttxuxxxututtt精确解为txcos二：数值解法思想和步骤2.1:网格剖分为了用差分方法求解上述问题，将求解区域(,)|01,01xtxt作剖分。将空间区间[0,1]作m等分，将时间[0,1]区间作n等分，并记1/,1/,,0,,0jkhmnxjhjmtkkn。分别称h和为空间和时间步长。用两簇平行直线,0,,0jkxxjmttkn将分割成矩形网格。2.2：差分格式的建立0uutx2.2.1：Lax-Friedrichs方法对时间、空间采用中心差分使得2h11111)(21uuxuuuuutukjkjkjkjkjkj则由上式得到Lax-Friedrichs格式111111()202kkkkkjjjjjuuuuuh截断误差为()[]kkkjhjjRuLuLu111111()22kkkkkkkjjjjjjjuuuuuuuhtx232223(),(0,0)26kkjjuuhOhjmkntx所以Lax-Friedrichs格式的截断误差的阶式2()Oh令/sh：则可得差分格式为111111(),(0,0)222kkkkkjjjjjssuuuuujmkn0cos()(0)jjuxjm0cos(),cos(),(0)kkkmkututkn其传播因子为：eeGhihishihiee221，化简可得：hsGhishGsin11,sincos,222所以当1s时，1,G,格式稳定。*2.2.2：LaxWendroff方法用牛顿二次插值公式可以得到LaxWendroff的差分格式，在此不详细分析，它的截断误差为h22，是二阶精度；当2s时，1,G，格式稳定。在这里主要用它与上面一阶精度的Lax-Friedrichs方法进行简单对比。2.3差分格式的求解因为1s时格式稳定，不妨取1/90h1/100，则s=0.9差分格式111111(),(0,0)222kkkkkjjjjjssuuuuujmkn写成如下矩阵形式：10111212111100000222211000002222110000002222000000000000000kkkkkkmmkmkmssuussuuuussuu则需要通过对k时间层进行矩阵作用求出k+1时间层。对上面的矩阵形式通过matlab编出如附录的程序求出数值解、真实解和误差。2.5算法以及结果function[PUExt]=PDEHyperbolic(uX,uT,M,N,C,type)formatlong%一阶双曲型方程的差分格式%[PUExt]=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)%方程:u_t+C*u_x=00=t=uT,0=x=uX%初值条件:u(x,0)=phi(x)%输出参数:U-解矩阵%x-横坐标%t-纵坐标，时间%输入参数:uX-变量x的上界%uT-变量t的上界%M-变量x的等分区间数%N-变量t的等分区间数%C-系数%phi-初值条件函数，定义为内联函数%psi1,psi2-边值条件函数，定义为内联函数%type-差分格式，从下列值中选取%-type='LaxFriedrichs'，采用Lax-Friedrichs差分格式求解%-type='LaxWendroff'，采用Lax-Wendroff差分格式求解h=uX/M;%变量x的步长k=uT/N;%变量t的步长r=k/h;%步长比x=(0:M)*h;t=(0:N)*k;U=zeros(M+1,N+1);%初值条件fori=1:M+1U(i,1)=cos(pi*x(i));P(i,1)=cos(pi*x(i));E(i,1)=0;end%边值条件forj=1:N+1U(1,j)=cos(pi*t(j));E(1,j)=0;P(1,j)=cos(pi*t(j));U(M+1,j)=-cos(pi*t(j));P(M+1,j)=-cos(pi*t(j));E(M+1,j)=0;endswitchtypecase'LaxFriedrichs'ifabs(C*r)1disp('|C*r|1，Lax-Friedrichs差分格式不稳定！')end%逐层求解forj=1:Nfori=2:MU(i,j+1)=(U(i+1,j)+U(i-1,j))/2-C*r*(U(i+1,j)-U(i-1,j))/2;P(i,j+1)=cos(pi*(x(i)+t(j+1)));E(i,j+1)=abs(U(i,j+1)-cos(pi*(x(i)+t(j+1))));endend%Lax-Wendroff差分格式case'LaxWendroff'ifabs(C*r)1disp('|C*r|1，Lax-Wendroff差分格式不稳定！')end%逐层求解forj=1:Nfori=2:MU(i,j+1)=U(i,j)-C*r*(U(i+1,j)-U(i-1,j))/2+C^2*r^2*(U(i+1,j)-2*U(i,j)+U(i-1,j))/2;P(i,j+1)=cos(pi*(x(i)+t(j+1)));E(i,j+1)=abs(U(i,j+1)-cos(pi*(x(i)+t(j+1))));endendotherwisedisp('差分格式类型输入有误！')return;endU=U';P=P';E=E';%作出图形精确解mesh(x,t,P);title('一阶双曲型方程的精确解图像');xlabel('空间变量x');ylabel('时间变量t');zlabel('一阶双曲型方程的解P')%作出图形数值解mesh(x,t,U);title([type'格式求解一阶双曲型方程的解的图像']);xlabel('空间变量x');ylabel('时间变量t');zlabel('一阶双曲型方程的解U')return;命令窗口输入：uX=1;uT=1;M=90;N=100;C=-1;phi=inline('cos(pi*x)');psi1=inline('cos(pi*t)');psi2=inline('-cos(pi*t)');type='LaxFriedrichs'或type='LaxWendroff';[PUExt]=PDEHyperbolic(uX,uT,M,N,C,type)从matlab的数值解法结果中抽出一部分数据进行比较表1LaxFriedrichs格式jk（x,t）数值解真实解误差4611(0.5,0.1)-0.308981-0.3090170.0000364621(0.5,0.2)-0.587647-0.5877850.0001384631(0.5,0.3)-0.808731-0.8090170.0002864641(0.5,0.4)-0.950609-0.9510560.0004484651(0.5,0.5)-0.999409-1.0000000.0005914661(0.5,0.6)-0.950496-0.9510570.0005604671(0.5,0.7)-0.808539-0.8090170.0004784681(0.5,0.8)-0.587437-0.5877850.0003484691(0.5,0.9)-0.308833-0.3090170.00018446101(0.5,1.0)-0.000002-0.0000000.000002表2LaxWendroff格式备注：本来100k090j0，，但是由于matlab中下标必须从大于0开始，所以在程序中101k191j1，图像分析:jk（x,t）数值解真实解误差4611(0.5,0.1)-0.309005-0.3090170.0000124621(0.5,0.2)-0.587765-0.5877850.0000204631(0.5,0.3)-0.808995-0.8090170.0000224641(0.5,0.4)-0.951040-0.9510560.0000164651(0.5,0.5)-0.999999-1.0000000.0000014661(0.5,0.6)-0.951074-0.9510570.0000174671(0.5,0.7)-0.809051-0.8090170.0000344681(0.5,0.8)-0.587833-0.5877850.0000484691(0.5,0.9)-0.309074-0.3090170.00005746101(0.5,1.0)-0.000006-0.0000000.000006结果分析：从表1和表2可以看出LaxFriedrichs格式和LaxWendroff格式的真值得误差都比较小，而LaxWendroff格式虽然精度比LaxFriedrichs的精度高，但是在网格点划分比较细的情况下，二者的差别不大。从三个图像的结果看出，二者都拟合的相当好，并且结果都稳定。