《偏微分方程数值解法MATLAB源码》由会员分享，可在线阅读，更多相关《偏微分方程数值解法MATLAB源码(27页珍藏版)》请在人人文库网上搜索。

1、源码【更新完毕】偏微分方程数值解法的MATLAB原创 说明：由于偏微分的程序都比较长，比其他的算法稍复杂一些，所以另开一贴，专门上传偏微分的程序 谢谢大家的支持！ 其他的数值算法见：./Announce/Announce.asp?BoardID=209&id=8245004 、古典显式格式求解抛物型偏微分方程(一维热传导方程)1 function U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) %古典显式格式求解抛物型偏微分方程 %U x t=PDEParabolicClassicalExplicit(uX,uT,p。

2、hi,psi1,psi2,M,N,C) % %方程：u_t=C*u_xx 0 if r ) 不稳定 0.5, disp(r end 计算初值和边值%U=zeros(M+1,N+1); i=1:M+1 for U(i,1)=phi(x(i); end j=1:N+1 for U(1,j)=psi1(t(j); U(M+1,j)=psi2(t(j); end 逐层求解%j=1:N for i=2:M for U(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j); end end U=U; %作出图形mesh(x,t,U); ) 古典显式格式，一维热传导方程的解的图像ti。

3、tle(x) xlabel(空间变量t) 时间变量 ylabel(U) zlabel(一维热传导方程的解 return; 古典显式格式不稳定情况2 / 16 古典显式格式稳定情况2、古典隐式格式求解抛物型偏微分方程(一维热传导方程) function U x t=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C) %古典隐式格式求解抛物型偏微分方程 %U x t=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C) % %方程：u_t=C*u_xx 0 1 if 13 / 。

4、16 ) 差分格式不稳定！，Lax-Friedrichs disp(|C*r|1 end 逐层求解 % j=1:N for i=2:M for U(i,j+1)=(U(i+1,j)+U(i-1,j)/2-C*r*(U(i+1,j)-U(i-1,j)/2; end end %Courant-Isaacson-Rees差分格式 CourantIsaacsonRees case C0 C*r1 if ) disp(Courant-Isaacson-Lees差分格式不稳定！ end %逐层求解 j=1:N for i=2:M for U(i,j+1)=C*r*U(i-1,j)+(1-C*r)*U(i,。

5、j); end end end %Leap-Frog(蛙跳)差分格式 LeapFrog case psi2=); 请输入第二层初值条件函数： phi2=input( abs(C*r)1 if ) Leap-Frog差分格式不稳定！ disp(|C*r|1， end %第二层初值条件 i=1:M+1 for U(i,2)=phi2(x(i); 14 / 16 end %逐层求解 j=2:N for i=2:M for U(i,j+1)=U(i,j-1)-C*r*(U(i+1,j)-U(i-1,j); end end 差分格式 %Lax-WendroffLaxWendroff case abs(C。

6、*r)1 if ) 差分格式不稳定！disp(|C*r|1，Lax-Wendroff end 逐层求解 % j=1:N for i=2:M for U(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; end end %Crank-Nicolson隐式差分格式，需调用追赶法求解三对角线性方程组的算法 CrankNicolson case Diag=zeros(1,M-1);%矩阵的对角线元素 Low=zeros(1,M-2);%矩阵的下对角线元素 Up=zeros(1,M-2);%矩阵的上对角。

7、线元素 i=1:M-2 for Diag(i)=4; Low(i)=-r*C; Up(i)=r*C; end Diag(M-1)=4; B=zeros(M-1,M-1); i=1:M-2 for B(i,i)=4; B(i,i+1)=-r*C; B(i+1,i)=r*C; end B(M-1,M-1)=4; ) %逐层求解，需要使用追赶法(调用函数EqtsForwardAndBackward j=1:N for b1=zeros(M-1,1); b1(1)=r*C*(U(1,j+1)+U(1,j)/2; 15 / 16 b1(M-1)=-r*C*(U(M+1,j+1)+U(M+1,j)/2; b=B*U(2:M,j)+b1; U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b); end otherwise ) 差分格式类型输入有误！ disp( return; end U=U; 作出图形%mesh(x,t,U); ); 格式求解一阶双曲型方程的解的图像title(type x); 空间变量 xlabel(t); ylabel(时间变量U); 一阶双曲型方程的解 zlabel( return; 16 / 16。