雙曲方程基于matlab的數(shù)值解法

./雙曲型方程基于MATLAB的數(shù)值解法〔數(shù)學(xué)1201,曉云,41262022一：一階雙曲型微分方程的初邊值問(wèn)題精確解為二：數(shù)值解法思想和步驟2.1:網(wǎng)格剖分為了用差分方法求解上述問(wèn)題,將求解區(qū)域作剖分。將空間區(qū)間作等分,將時(shí)間區(qū)間作等分,并記。分別稱(chēng)和為空間和時(shí)間步長(zhǎng)。用兩簇平行直線(xiàn)將分割成矩形網(wǎng)格。2.2：差分格式的建立2.2.1：Lax-Friedrichs方法對(duì)時(shí)間、空間采用中心差分使得則由上式得到Lax-Friedrichs格式截?cái)嗾`差為所以Lax-Friedrichs格式的截?cái)嗾`差的階式令：則可得差分格式為其傳播因子為：化簡(jiǎn)可得：所以當(dāng)時(shí),,格式穩(wěn)定。*2.2.2：LaxWendroff方法用牛頓二次插值公式可以得到LaxWendroff的差分格式,在此不詳細(xì)分析,它的截?cái)嗾`差為,是二階精度；當(dāng)時(shí),,格式穩(wěn)定。在這里主要用它與上面一階精度的Lax-Friedrichs方法進(jìn)行簡(jiǎn)單對(duì)比。2.3差分格式的求解因?yàn)闀r(shí)格式穩(wěn)定,不妨取,則s=0.9差分格式寫(xiě)成如下矩陣形式：則需要通過(guò)對(duì)k時(shí)間層進(jìn)行矩陣作用求出k+1時(shí)間層。對(duì)上面的矩陣形式通過(guò)matlab編出如附錄的程序求出數(shù)值解、真實(shí)解和誤差。2.5算法以及結(jié)果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>%輸出參數(shù):U-解矩陣%x-橫坐標(biāo)%t-縱坐標(biāo),時(shí)間%輸入?yún)?shù):uX-變量x的上界%uT-變量t的上界%M-變量x的等分區(qū)間數(shù)%N-變量t的等分區(qū)間數(shù)%C-系數(shù)%phi-初值條件函數(shù),定義為聯(lián)函數(shù)%psi1,psi2-邊值條件函數(shù),定義為聯(lián)函數(shù)%type-差分格式,從下列值中選取%-type='LaxFriedrichs',采用Lax-Friedrichs差分格式求解%-type='LaxWendroff',采用Lax-Wendroff差分格式求解h=uX/M;%變量x的步長(zhǎng)k=uT/N;%變量t的步長(zhǎng)r=k/h;%步長(zhǎng)比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差分格式不穩(wěn)定！'>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差分格式不穩(wěn)定！'>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<'差分格式類(lèi)型輸入有誤！'>return;endU=U';P=P';E=E';%作出圖形精確解mesh<x,t,P>;title<'一階雙曲型方程的精確解圖像'>;xlabel<'空間變量x'>;ylabel<'時(shí)間變量t'>;zlabel<'一階雙曲型方程的解P'>%作出圖形數(shù)值解mesh<x,t,U>;title<[type'格式求解一階雙曲型方程的解的圖像']>;xlabel<'空間變量x'>;ylabel<'時(shí)間變量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的數(shù)值解法結(jié)果中抽出一部分?jǐn)?shù)據(jù)進(jìn)行比較表1LaxFriedrichs格式j(luò)k〔x,t數(shù)值解真實(shí)解誤差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格式j(luò)k〔x,t數(shù)值解真實(shí)解誤差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備注：本來(lái)