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

1、双曲型方程基于 MATLAB 的数值解法（数学 1201, 晓云 ,41262022）：一阶双曲型微分方程的初边值问题uu0,0 x 1,0 t 1.txu( x,0) cos( x),0 x 1. 精确解为 cos x t u(0, t) u(1,t) cos( t),0 t 1.二：数值解法思想和步骤2.1: 网格剖分为了用差分方法求解上述问题，将求解区域(x,t ) |0 x 1,0 t 1作剖分。将空间区间 0,1 作 m 等分，将时间0,1 区间作 n 等分，并记h 1/ m,1/ n, xjjh ,0jm, tk k ,0 kn。分别称 h 和 为空间和时间步长。用两簇平行直线x

2、xj ,0j m, ttk,0k n 将 分割成矩形网格。2.2 ：差分格式的建立u0xut2.2.1 ： Lax-Friedrichs 方法对时间、空间采用中心差分使得k1 ujkujk112 (u jku j 1)kku j 1 u j 12hux则由上式得到 Lax-Friedrichs 格式k 1 1 k k uj 2 (uj 1 uj 1)kku kj 1 ukj 1 02h截断误差为Rjk (u) Lhukj Lukjkukj1k2(uj 1ujk 1)k k k kj 1ukj 1 ukj 1ujkukj2h t x2 k 2ujh23kujO( h2),(0 j m,0 k n

3、)2t2 63 x所以 Lax-Friedrichs 格式的截断误差的阶式 O( h2)令 s /h ：则可得差分格式为k1ukj 1s k 1 k2uj 1 2(uj 1k s k ujk 1) 2ukj1,(0m,0n)0u0jcos( xj )(0 jm)ku0cos( tk ),umkcos(tk ),(0k n)其传播因子为：所以当 sihecos2111 时， G* 2.2.2 ： LaxWendroffsi2eissin h2sin h1, 格式稳定方法用牛顿二次插值公式可以得到LaxWendroff化简可得：的差分格式，在此不详细分析，它的截断误差为h ，是二阶精度；当 s 2

4、 时，G ,1，格式稳定。在这里主要用它与上面一阶精度的 Lax-Friedrichs 方法进行简单对比2.3 差分格式的求解因为 s 1时格式稳定，不妨取 h 11/19000 ，则 s=0.9差分格式 ukj 1 sukj 1 1(ukj 1 ukj 1) sukj 1,(0 j m,0 k n)写成如下矩阵形式：1s01s0222201s01s2222MMMM00000000L000ku0L000ku1MMMMM k umk 2L1s01suk um 12222L000umku1u2k1m1000L000则需要通过对k 时间层进行矩阵作用求出k+1 时间层对上面的矩阵形式通过 matla

5、b 编出如附录的程序求出数值解、真实解和误差2.5 算法以及结果function P U E x t=PDEHyperbolic(uX,uT,M,N,C,type)format long%一阶双曲型方程的差分格式%P U E x t=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)% 方程 :u_t+C*u_x=00 = t = uT, 0 = x 1disp( |C*r|1 ， Lax-Friedrichs 差分格式不稳定！ )end% 逐层求解for j=1:Nfor i=2:MU(i,j+1)=(U(i+1,j)+U(i-1,j)/2-C*r*(

6、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 LaxWendroffif abs(C*r)1disp( |C*r|1 ， Lax-Wendroff 差分格式不稳定！ )end% 逐层求解for j=1:Nfor i=2:MU(i,j+1)=U(i,j)-C*r*(U(i+1,j)-U(i-1,j)/2+C2*r2*(U(i+1,j)-2*U(i,j)+U(i-1,j)/2;P(i,j+1)=cos(

7、pi*(x(i)+t( j+1);j+1);E(i,j+1)=abs(U(i,j+1)-cos(pi*(x(i)+t(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 ）

8、; 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(-c os(pi*t);type=LaxFriedrichs 或 type=LaxWendroff;P U E x t=PDEHyperbolic(uX,uT,M,N,C,type)从 matlab 的数值解法结果中抽出一部分数据进行比较表1LaxFriedrichs格式jk（x,t）数值解真实解误差4611(0.5,0.1)-0.308981-0.3090170

9、.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.3088

10、33-0.3090170.00018446101(0.5,1.0)-0.000002-0.0000000.000002表2LaxWendroff格式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.

11、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备注：本来 0 j 90，0 k 100 ，但是由于 matlab 中下标必须从大于 0 开始，所以在程序中 1 j 91，1 k 101图像分析:Figure 1文件（D 糯 （Y）插入Q）二具CD 黨面9）京口包D 弟助凹曰巳齡|r巴働凰彳彳尿|口eTTq-一阶双曲型方程的精确解團像时间变量t0 0空间变量X一阶强曲型万程BWUHFigue 1 m苗噩6sg 多eH0 細旦 團口国确蓉巨 口鸟Igllrl-P-R8lellllftllxld-IPI-1IEII-LaxFriedrichs 诙n#另记HHfltlDsj?