![偏微分方程数值解法的MATLAB源码参考模板_第1页](http://file3.renrendoc.com/fileroot_temp3/2022-3/5/7313e52b-a95c-42da-a82e-b8bb2fb123ea/7313e52b-a95c-42da-a82e-b8bb2fb123ea1.gif)
![偏微分方程数值解法的MATLAB源码参考模板_第2页](http://file3.renrendoc.com/fileroot_temp3/2022-3/5/7313e52b-a95c-42da-a82e-b8bb2fb123ea/7313e52b-a95c-42da-a82e-b8bb2fb123ea2.gif)
![偏微分方程数值解法的MATLAB源码参考模板_第3页](http://file3.renrendoc.com/fileroot_temp3/2022-3/5/7313e52b-a95c-42da-a82e-b8bb2fb123ea/7313e52b-a95c-42da-a82e-b8bb2fb123ea3.gif)
![偏微分方程数值解法的MATLAB源码参考模板_第4页](http://file3.renrendoc.com/fileroot_temp3/2022-3/5/7313e52b-a95c-42da-a82e-b8bb2fb123ea/7313e52b-a95c-42da-a82e-b8bb2fb123ea4.gif)
![偏微分方程数值解法的MATLAB源码参考模板_第5页](http://file3.renrendoc.com/fileroot_temp3/2022-3/5/7313e52b-a95c-42da-a82e-b8bb2fb123ea/7313e52b-a95c-42da-a82e-b8bb2fb123ea5.gif)
版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
1、原创偏微分方程数值解法的MATLAB源码【更新完毕】说明:由于偏微分的程序都比较长,比其他的算法稍复杂一些,所以另开一贴,专门上传偏微分的程序谢谢大家的支持!其他的数值算法见:./Announce/Announce.asp?BoardID=209&id=82450041、古典显式格式求解抛物型偏微分方程(一维热传导方程)1 / 23function U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%古典显式格式求解抛物型偏微分方程%U x t=PDEParabolicClassicalExplicit(uX,u
2、T,phi,psi1,psi2,M,N,C)%方程:u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT%初值条件:u(x,0)=phi(x)%边值条件:u(0,t)=psi1(t), u(uX,t)=psi2(t)%输出参数:U -解矩阵,第一行表示初值,第一列和最后一列表示边值,第二行表示第2层% x -空间变量% t -时间变量%输入参数:uX -空间变量x的取值上限%
3、160; uT -时间变量t的取值上限% phi -初值条件,定义为内联函数% psi1 -边值条件,定义为内联函数% psi2 -边值条件,定义为内联函数% M -沿x轴的等分区间数% N -沿t轴的等分区间数% &
4、#160; C -系数,默认情况下C=1%应用举例:%uX=1;uT=0.2;M=15;N=100;C=1;%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');%U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C);%设置参数C的默认值if nargin=7 C=1;end%计算步长dx=uX/M;%x的步长dt=uT/N;%t的步长x=(0:M)*dx;t=(0:N)*d
5、t;r=C*dt/dx/dx;%步长比r1=1-2*r;if r > 0.5 disp('r > 0.5,不稳定')end%计算初值和边值U=zeros(M+1,N+1);for i=1:M+1 U(i,1)=phi(x(i);endfor j=1:N+1 U(1,j)=psi1(t(j); U(M+1,j)=psi2(t(j);end%逐层求解for j=1:N for i=2:M &
6、#160;U(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j); endendU=U'%作出图形mesh(x,t,U);title('古典显式格式,一维热传导方程的解的图像')xlabel('空间变量 x')ylabel('时间变量 t')zlabel('一维热传导方程的解 U')return;古典显式格式不稳定情况古典显式格式稳定情况2、古典隐式格式求解抛物型偏微分方程(一维热传导方程)function U x t=PDEParabolicClassicalImpli
7、cit(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 <= x <= uX,0 <= t <= uT%初值条件:u(x,0)=phi(x)%边值条件:u(0,t)=psi1(t), u(uX,t)=psi2(t)%输出参数:U -解矩阵,第一行表示初值,第一列和最后一列表示边值,第二行表示第2层% x -空间
8、变量% t -时间变量%输入参数:uX -空间变量x的取值上限% uT -时间变量t的取值上限% phi -初值条件,定义为内联函数% psi1 -边值条件,定义为内联函数% psi2 -边值条件,定义为内联函数% &
9、#160; M -沿x轴的等分区间数% N -沿t轴的等分区间数% C -系数,默认情况下C=1%应用举例:%uX=1;uT=0.2;M=50;N=50;C=1;%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');%U x t=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C);%设置参数C
10、的默认值if nargin=7 C=1;end%计算步长dx=uX/M;%x的步长dt=uT/N;%t的步长x=(0:M)*dx;t=(0:N)*dt;r=C*dt/dx/dx;%步长比Diag=zeros(1,M-1);%矩阵的对角线元素Low=zeros(1,M-2);%矩阵的下对角线元素Up=zeros(1,M-2);%矩阵的上对角线元素for i=1:M-2 Diag(i)=1+2*r; Low(i)=-r; Up(i)=-r;endDiag(M-1)=1+2*r;%计算初值和边值U
11、=zeros(M+1,N+1);for i=1:M+1 U(i,1)=phi(x(i);endfor j=1:N+1 U(1,j)=psi1(t(j); U(M+1,j)=psi2(t(j);end%逐层求解,需要使用追赶法(调用函数EqtsForwardAndBackward)for j=1:N b1=zeros(M-1,1); b1(1)=r*U(1,j+1); b1(M-1)=r*U(M+1,j+1); b=
12、U(2:M,j)+b1; U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);endU=U'%作出图形mesh(x,t,U);title('古典隐式格式,一维热传导方程的解的图像')xlabel('空间变量 x')ylabel('时间变量 t')zlabel('一维热传导方程的解 U')return;此算法需要使用追赶法求解三对角线性方程组,这个算法在上一篇帖子中已经给出,为了方便,再给出来追赶法解三对角线性方程组function x=EqtsForw
13、ardAndBackward(L,D,U,b)%追赶法求解三对角线性方程组Ax=b%x=EqtsForwardAndBackward(L,D,U,b)%x:三对角线性方程组的解%L:三对角矩阵的下对角线,行向量%D:三对角矩阵的对角线,行向量%U:三对角矩阵的上对角线,行向量%b:线性方程组Ax=b中的b,列向量%应用举例:%L=-1 -2 -3;D=2 3 4 5;U=-1 -2 -3;b=6 1 -2 1'%x=EqtsForwardAndBackward(L,D,U,b)%检查参数的输入是否正确n=length(D);m=length(b);n1=length(L);n2=len
14、gth(U);if n-n1 = 1 | n-n2 = 1 | n = m disp('输入参数有误!') x=' ' return;end%追的过程for i=2:n L(i-1)=L(i-1)/D(i-1); D(i)=D(i)-L(i-1)*U(i-1);endx=zeros(n,1);x(1)=b(1);for i=2:n x(i)=b(i)-L(i-1)*x(i-1);end%赶的过程x(n)=x(
15、n)/D(n);for i=n-1:-1:1 x(i)=(x(i)-U(i)*x(i+1)/D(i);endreturn;古典隐式格式在以后的程序中,我们都取C=1,不再作为一个输入参数处理3、Crank-Nicolson隐式格式求解抛物型偏微分方程需要调用追赶法的程序function U x t=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)%Crank-Nicolson隐式格式求解抛物型偏微分方程%U x t=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)%方程:u_t=u_xx 0 <= x
16、<= uX,0 <= t <= uT%初值条件:u(x,0)=phi(x)%边值条件:u(0,t)=psi1(t), u(uX,t)=psi2(t)%输出参数:U -解矩阵,第一行表示初值,第一列和最后一列表示边值,第二行表示第2层% x -空间变量% t -时间变量%输入参数:uX -空间变量x的取值上限% uT -时间变量t的取值上限%
17、 phi -初值条件,定义为内联函数% psi1 -边值条件,定义为内联函数% psi2 -边值条件,定义为内联函数% M -沿x轴的等分区间数% N -沿t轴的等分区间数%应用举例:%uX=1;uT=0.2;M=50;N=50;%phi=inline('sin(pi*x)');psi1=inl
18、ine('0');psi2=inline('0');%U x t=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N);%计算步长dx=uX/M;%x的步长dt=uT/N;%t的步长x=(0:M)*dx;t=(0:N)*dt;r=dt/dx/dx;%步长比Diag=zeros(1,M-1);%矩阵的对角线元素Low=zeros(1,M-2);%矩阵的下对角线元素Up=zeros(1,M-2);%矩阵的上对角线元素for i=1:M-2 Diag(i)=1+r; Low(i)=-r/2;
19、 Up(i)=-r/2;endDiag(M-1)=1+r;%计算初值和边值U=zeros(M+1,N+1);for i=1:M+1 U(i,1)=phi(x(i);endfor j=1:N+1 U(1,j)=psi1(t(j); U(M+1,j)=psi2(t(j);endB=zeros(M-1,M-1);for i=1:M-2 B(i,i)=1-r; B(i,i+1)=r/2; B(i+1,i)=r/2;endB
20、(M-1,M-1)=1-r;%逐层求解,需要使用追赶法(调用函数EqtsForwardAndBackward)for j=1:N b1=zeros(M-1,1); b1(1)=r*(U(1,j+1)+U(1,j)/2; b1(M-1)=r*(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);endU=U'%作出图形mesh(x,t,
21、U);title('Crank-Nicolson隐式格式,一维热传导方程的解的图像')xlabel('空间变量 x')ylabel('时间变量 t')zlabel('一维热传导方程的解 U')return;Crank-Nicolson隐式格式4、正方形区域Laplace方程Diriclet问题的求解需要调用Jacobi迭代法和Guass-Seidel迭代法求解线性方程组function U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)%正方形
22、区域Laplace方程的Diriclet边值问题的差分求解%此程序需要调用Jacobi迭代法或者Guass-Seidel迭代法求解线性方程组%U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)%方程:u_xx+u_yy=0 0<=x,y<=ub%边值条件:u(0,y)=phi1(y)% u(ub,y)=phi2(y)% u(x,0)=psi1(x)% &
23、#160; u(x,ub)=psi2(x)%输出参数:U -解矩阵,第一行表示y=0时的值,第二行表示第y=h时的值% x -横坐标% y -纵坐标%输入参数:ub -变量边界值的上限% phi1,phi2,psi1,psi2 -边界函数,定义为内联函数% M -横纵坐标的等分区间数
24、% type -求解差分方程的迭代格式,若type='Jacobi',采用Jacobi迭代格式% 若type='GS',采用Guass-Seidel迭代格式。默认情况下,type='GS'%应用举例:%ub=4;M=20;%phi1=inline('y*(4-y)');phi2=inline('0');psi1=inline('
25、;sin(pi*x/4)');psi2=inline('0');%U x y=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,'GS');if nargin=6 type='GS'end%步长h=ub/M;%横纵坐标x=(0:M)*h;y=(0:M)*h;%差分格式的矩阵形式AU=K%构造矩阵AM2=(M-1)2;A=zeros(M2);for i=1:M2 A(i,i)=4;endfor i=1:M2-1
26、 if mod(i,M-1)=0 A(i,i+1)=-1; A(i+1,i)=-1; endendfor i=1:M2-M+1 A(i,i+M-1)=-1; A(i+M-1,i)=-1;endU=zeros(M+1);%边值条件for i=1:M+1 U(i,1)=psi1(i-1)*h); U(i,M+1)=psi2(i-1)
27、*h); U(1,i)=phi1(i-1)*h); U(M+1,i)=phi2(i-1)*h);end%构造KK=zeros(M2,1);for i=1:M-1 K(i)=U(i+1,1); K(M2-i+1)=U(i+1,M+1);endK(1)=K(1)+U(1,2);K(M-1)=K(M-1)+U(M+1,2);K(M2-M+2)=K(M2-M+2)+U(1,M);K(M2)=K(M2)+U(M+1,M);for i=2:M-2 K(M-1)*(i-1)+1)
28、=U(1,i+1); K(M-1)*i)=U(M+1,i+1);endx0=ones(M2,1);switch type %调用Guass-Seidel迭代法求解线性方程组AU=K case 'Jacobi' X=EqtsJacobi(A,K,x0); %调用Guass-Seidel迭代法求解线性方程组AU=K case 'GS' &
29、#160; X=EqtsGS(A,K,x0); otherwise disp('差分格式类型输入错误') return;end%把求解结果化成矩阵型式for i=2:M for j=2:M U(j,i)=X(j-1+(M-1)*(i-2); endendU=U'%作出图形mesh(x,y,U)
30、;title('五点差分格式Laplace方程Diriclet问题的解的图像')xlabel('x')ylabel('y')zlabel('Laplace方程Diriclet问题的解 U')return;正方形区域Laplace方程五点差分格式5、一阶双曲型方程的差分方法function U x t=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)%一阶双曲型方程的差分格式%U x t=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)%方程:u_t
31、+C*u_x=0 0 <= t <= uT, 0 <= x <= uX%初值条件:u(x,0)=phi(x)%输出参数:U -解矩阵,第一行表示初值,第二行表示第2个时间层% x -横坐标% t -纵坐标,时间%输入参数:uX -变量x的上界% uT -变量t的上界% M -变量x的等分区间数%
32、; N -变量t的等分区间数% C -系数% phi -初值条件函数,定义为内联函数% psi1,psi2 -边值条件函数,定义为内联函数% type -差分格式,从下列值中选取% -type='LaxFriedrichs',采用Lax-Friedrichs差分格
33、式求解% -type='CourantIsaacsonRees',采用Courant-Isaacson-Rees差分格式求解% -type='LeapFrog',采用Leap-Frog(蛙跳)差分格式求解% -type='LaxWendroff
34、39;,采用Lax-Wendroff差分格式求解% -type='CrankNicolson',采用Crank-Nicolson差分格式求解,此格式需调用追赶法% 求解三对角线性方程组%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);%初值条件for i=1:M+1
35、60; U(i,1)=phi(x(i);end%边值条件for j=1:N+1 U(1,j)=psi1(t(j); U(M+1,j)=psi2(t(j); %U(1,j)=NaN; %U(M+1,j)=NaN;endswitch type %Lax-Friedrichs差分格式 case 'LaxFriedrichs' if abs(C*r)>
36、;1 disp('|C*r|>1,Lax-Friedrichs差分格式不稳定!') end %逐层求解 for j=1:N for i=2:M
37、60; 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差分格式 case 'CourantIsaacsonRees' &
38、#160; if C<0 disp('C<0,采用前差公式') if C*r<-1 disp('Courant-Isaacson-Lees差分格式不稳定!')
39、 end %逐层求解 for j=1:N for i=2:M U(i,j+1)=(1+C*r
40、)*U(i,j)-C*r*U(i+1,j); end end else disp('C>0,采用后差公式')
41、 if C*r>1 disp('Courant-Isaacson-Lees差分格式不稳定!') end %逐层求解 for j=1:N
42、; for i=2:M U(i,j+1)=C*r*U(i-1,j)+(1-C*r)*U(i,j); end end
43、60; end %Leap-Frog(蛙跳)差分格式 case 'LeapFrog' phi2=input('请输入第二层初值条件函数:psi2='); if abs(C*r)>1
44、160; disp('|C*r|>1,Leap-Frog差分格式不稳定!') end %第二层初值条件 for i=1:M+1 U(i,2)=phi2(x(i); end
45、60;%逐层求解 for j=2:N for i=2:M U(i,j+1)=U(i,j-1)-C*r*(U(i+1,j)-U(i-1,j); end end&
46、#160; %Lax-Wendroff差分格式 case 'LaxWendroff' if abs(C*r)>1 disp('|C*r|>1,Lax-Wendroff差分格式不稳定!') end
47、0; %逐层求解 for j=1:N for i=2:M 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隐式差分格式,需调用追赶法求解三对角线性方程组的算法 case 'CrankNicolson' Diag=zeros(1,M-1);%矩阵的对角线元素
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 新版北师大买文具教学设计心得体会
- 常见燃料的优缺点分析
- 苏教版五年级生字表教学设计教学方法
- 会计服务合同模板
- 进口三方协议的范本(标准版)
- 多人合伙退伙协议(标准版)
- 八年级数学人教版(上册)期末测试(A卷)
- 八年级语文人教部编版(上册)第5单元 总结复习
- 六年级语文人教部编版(上册)8. 灯光
- 统编版四年级语文上册第3课《现代诗二首》精美课件
- EV71疫苗质量控制和评价技术要点
- GB/T 7036.2-2017充气轮胎内胎第2部分:摩托车轮胎内胎
- GB/T 38530-2020城镇液化天然气(LNG)气化供气装置
- 小学英语外研版(三起点)五年级下册全册课文翻译(1-10模块)
- 高中化学竞赛有机化学基础课件
- 浙江省名校协作体2022-2023学年高二上学期开学考试数学试卷及答案
- 《中医筋伤学教材》word版参考模板
- 海尔xqg90-bdx1226滚筒洗衣机中文服务手册
- 初中劳动人教七年级上册目录 劳技课第一课《日常礼仪》PPT
- 园林植物育种学-03种质资源
- A修标准项目专业资料
评论
0/150
提交评论