偏微分方程数值解法的MATLAB源码.docx
![资源得分’ title=](/images/score_1.gif)
![资源得分’ title=](/images/score_1.gif)
![资源得分’ title=](/images/score_1.gif)
![资源得分’ title=](/images/score_1.gif)
![资源得分’ title=](/images/score_05.gif)
《偏微分方程数值解法的MATLAB源码.docx》由会员分享,可在线阅读,更多相关《偏微分方程数值解法的MATLAB源码.docx(6页珍藏版)》请在淘文阁 - 分享文档赚钱的网站上搜索。
1、精选优质文档-倾情为你奉上原创偏微分方程数值解法的MATLAB源码【更新完毕】说明:由于偏微分的程序都比较长,比其他的算法稍复杂一些,所以另开一贴,专门上传偏微分的程序谢谢大家的支持!其他的数值算法见:1、古典显式格式求解抛物型偏微分方程(一维热传导方程)function U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%古典显式格式求解抛物型偏微分方程%U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%方程:u_t=C*u_xx 0 = x
2、= uX,0 = t 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 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(一维热传
3、导方程的解 U)return;古典显式格式不稳定情况古典显式格式稳定情况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 = x = uX,0 = t = uT%初值条件:u(x,0)=phi(x)%边值条件:u(0,t)=psi1(t), u(uX,t)=p
4、si2(t)%输出参数:U -解矩阵,第一行表示初值,第一列和最后一列表示边值,第二行表示第2层% x -空间变量% t -时间变量%输入参数:uX -空间变量x的取值上限% uT -时间变量t的取值上限% phi -初值条件,定义为内联函数% psi1 -边值条件,定义为内联函数% psi2 -边值条件,定义为内联函数% 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=PDE
5、ParabolicClassicalImplicit(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)*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;%
6、计算初值和边值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%逐层求解,需要使用追赶法(调用函数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=U(2:M,j)+b1; U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);endU=U;%作出图形mesh(x,t,
7、U);title(古典隐式格式,一维热传导方程的解的图像)xlabel(空间变量 x)ylabel(时间变量 t)zlabel(一维热传导方程的解 U)return;此算法需要使用追赶法求解三对角线性方程组,这个算法在上一篇帖子中已经给出,为了方便,再给出来追赶法解三对角线性方程组function x=EqtsForwardAndBackward(L,D,U,b)%追赶法求解三对角线性方程组Ax=b%x=EqtsForwardAndBackward(L,D,U,b)%x:三对角线性方程组的解%L:三对角矩阵的下对角线,行向量%D:三对角矩阵的对角线,行向量%U:三对角矩阵的上对角线,行向量%b
8、:线性方程组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=length(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);
9、x(1)=b(1);for i=2:n x(i)=b(i)-L(i-1)*x(i-1);end%赶的过程x(n)=x(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(u
10、X,uT,phi,psi1,psi2,M,N)%方程:u_t=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的取值上限% uT -时间变量t的取值上限% phi -初值条件,定义为内联函数% psi1 -边值条件,定义为内联函数% psi2 -边值条件,定义为内联函数% M -沿x轴的等分区间数% N -沿t轴的等分区间数%应用举
11、例:%uX=1;uT=0.2;M=50;N=50;%phi=inline(sin(pi*x);psi1=inline(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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 解法 MATLAB 源码
![提示](https://www.taowenge.com/images/bang_tan.gif)
限制150内