《偏微分方程数值解法的MATLAB源码》由会员分享，可在线阅读，更多相关《偏微分方程数值解法的MATLAB源码（10页珍藏版）》请在金锄头文库上搜索。
1、原创偏微分方程数值解法的AAB源码【更新完毕】阐明:由于偏微分的程序都比较长，比其她的算法稍复杂某些，因此另开一贴,专门上传偏微分的程序谢谢人们的支持！其她的数值算法见：./Annce/Anounceap?BoI=209&id=245041、古典显式格式求解抛物型偏微分方程(一维热传导方程)function x t=PDEParboliClassilExplicit(uX,u,phi,psi1,ps2,M,C）古典显式格式求解抛物型偏微分方程U x t=DEarabolcClsscalpiit（uX,uT,phi，psi1,psi2,N,C）%方程:_=C*u_xx 0=x = uX,0 =
2、t = uT初值条件：(,)=ph（x)边值条件:u(,t)=psi1(t), u(uX，)=i2(t)%输出参数:U -解矩阵，第一行表达初值，第一列和最后一列表达边值,第二行表达第层 x 空间变量% t-时间变量%输入参数：X -空间变量x的取值上限% uT -时间变量的取值上限% hi-初值条件,定义为内联函数% p -边值条件,定义为内联函数% psi 边值条件,定义为内联函数 M -沿x轴的等分区间数 N -沿t轴的等分区间数% C -系数，默认状况下C=%应用举例:%uX;uT0.2;M=15;N=100;C1；pi=inne(sn(p)）;si1=inlne(0)；psi2=in
3、ne(0)；%U t=PParabolicClassaExpicit（uX,uT,hi,psi1，s2,M,N，C）;%设立参数C的默认值if nargi=7 =；nd计算步长dxX/M;%x的步长dtuN;t的步长x=(:M)dx；t=(0:N）*dt；r=Cdt/dxx;%步长比12*r;i r 0.disp（ 0.5，不稳定)end%计算初值和边值=zeros(M+1,N+1);ori=1:M+1U（i，1)=phi(i);edforj=1:+1U(1,）=ps1(j); U(M+1,j)=ps（t(j);nd%逐级求解r j=1：N fr=2：M U(i，j+)rU(i-1，j)+r1
4、*(i，）r*(i+,j); endend=U;%作出图形mesh(,t,U);title(古典显式格式，一维热传导方程的解的图像)labe（空间变量 x)ylbel(时间变量 t)zabel（一维热传导方程的解U)reurn;古典显式格式不稳定状况古典显式格式稳定状况2、古典隐式格式求解抛物型偏微分方程（一维热传导方程)funcinUx t=PDEPabolClasscalImlii（X,uT,p,psi1,p，M，N，C)%古典隐式格式求解抛物型偏微分方程U x t=DEParaboliCassicalImpici(u，u,ph,si1,i2，M,N,C)%方程:u_tC*ux= = uX
5、,0 = t= uT%初值条件:(,0）p(x)边值条件:u(,t)=psi1(t), u(,t)ps2(t)%输出参数：U -解矩阵，第一行表达初值，第一列和最后一列表达边值,第二行表达第2层% -空间变量% t -时间变量%输入参数:u-空间变量x的取值上限% -时间变量的取值上限% hi初值条件,定义为内联函数 pi1边值条件，定义为内联函数% ps2 边值条件,定义为内联函数% M-沿x轴的等分区间数 沿t轴的等分区间数% C-系数，默认状况下1%应用举例:%u;uT=.2；=50;N=0;C=1;%ph=inlne（sn(pi*);p1=inlie(0);si2=l(0)；% xt=
6、PDEaraolcClasicalIplit（uX，u,pi,psi1，psi,N,C)；%设立参数C的默认值f argin= C1;end计算步长dx=uXM;%x的步长tuTN;%的步长x(:M)*;t=（:N)*d；r=Cdt/dxdx；步长比Daes（1,M-1)；%矩阵的对角线元素Lowzeo（,M-2);%矩阵的下对角线元素=zers(1,M-2）;%矩阵的上对角线元素or i=1：M- ig(i)=1+*r; ow（i)-r; Up()=；endDiag(M)=*r;%计算初值和边值U=ers(M+1,N+);or=1：M （,1)=pi（x(i);enfor j=:N+1（1,
7、)(j）); U(1，j)=si(t(j)；ed逐级求解，需要使用追赶法(调用函数EqsForwadndakwar)for j=1:N b1=zers(1,); (1)=r*U(1,j+1)； b（M-1)=U(+1,j1); b=（2：M,j)b1； U（2:M,j+1)qsFrwadAndBcward(Lw，ag,U,);end=U;作出图形mesh（x，,U);titl(古典隐式格式,一维热传导方程的解的图像)xlabel(空间变量 x)ylabel(时间变量 t)zabel(一维热传导方程的解U)retun;此算法需要使用追赶法求解三对角线性方程组，这个算法在上一篇帖子中已经给出，为了
8、以便，再给出来追赶法解三对角线性方程组funon =EtsForrdAdacwad（L,D，U,b)%追赶法求解三对角线性方程组Ax%x=EtsrwdAndBcwad(L，D,U,b)%x:三对角线性方程组的解%L:三对角矩阵的下对角线,行向量%D：三对角矩阵的对角线,行向量%:三对角矩阵的上对角线,行向量%b:线性方程组Ax=b中的b，列向量%应用举例：%=- -2 -3;D=2 5;U=1 2 -3；b=6 -21;%x=EtFowarAndBackard（L,D,b)%检查参数的输入与否对的n=egth(D)；m=egt（b);1=len(L）;n2lengt（U）；f n-n =1 |
9、 n-2 1 | = m dsp(输入参数有误!) x ; return;n%追的过程for =2:n L(i-1)（-）/D(i-1); D()=D(i)-(i-1)U(-);endx=eros（n，1）;x()b();fo i=:n x(i)=(i)-L（i-1）*x(i-1）;end%赶的过程x(）x(n)D（n);fri=-:-1: x(i)（(i）-U（)*x(+1)D（i);endretur;古典隐式格式在后来的程序中，我们都取C=1,不再作为一种输入参数解决3、CraNclson隐式格式求解抛物型偏微分方程需要调用追赶法的程序unctionUx =PDEParabolicC（uX
10、,uT，p,psi1,si,M,N)Craicolso隐式格式求解抛物型偏微分方程U x tPDParbocN（X,h,psi1,pi,M,)方程：u_t=ux 0 = x uX,= t= uT%初值条件：u(x,0)phi()%边值条件：u(0,t）ps1(t), u(uX,t）psi2（t)%输出参数：U 解矩阵,第一行表达初值,第一列和最后一列表达边值，第二行表达第2层% x -空间变量% t -时间变量输入参数:X -空间变量x的取值上限 uT -时间变量t的取值上限% phi -初值条件,定义为内联函数% pi -边值条件,定义为内联函数% ps2 -边值条件，定义为内联函数% M
11、-沿x轴的等分区间数% N沿t轴的等分区间数%应用举例：uX=1;uT=0.2;5；=50;%hi=ili(si(pi*x)；psi1=inli(0）；pi2=inline(0);%U x t=DPaabolcCN(uX，u,phi,si,psi2，,N)；%计算步长duX/M;%x的步长duT；%t的步长x(:M）*dx;t(:)*t;=dtdx/dx;步长比iag=zeos(1,M1);矩阵的对角线元素ower（1,M-2）；%矩阵的下对角线元素Ueros(1,M-);%矩阵的上对角线元素for i=1:M-2 Diag(i)=+; Lo(i)=r/2;Up（i）-r/2;eD(M-1)=1+r；%计算初值和边值Uzeros（1,N+1）;ori=1:M+1 U（i,)=phi(i）)；nf j=1：N+1 U(1,j）=psi(t(j）)； U（M+1,j)=2(t(j)）;endBzes(-，-)；fo i=1:M-2 B(,i）1-r; （,+1）r2; B(+1,i)=r/2;endB(1，M-)=-r;%逐级求解,需要使用追赶法(调用函数EqtsFrardndBackwr
