最新偏微分方程的matlab解法
偏微分方程Matlab数值解法(补充4)
偏微分方程Matlab 数值解法（补充4）Matlab 可以求解一般的偏微分方程，也可以利用偏微分方程工具箱中给出的函数求解一些偏微分方程。
1 偏微分方程组求解Matlab 语言提供了pdepe()函数，可以直接求解偏微分方程(,,,)[(,,,)](,,,)m mu u u u C x t u x x f x t u s x t u x t x x x-∂∂∂∂∂=+∂∂∂∂∂ （4.1）这样，偏微分方程可以编写为以下函数的描述，其入口为[,,](,,,)x c f s pdefun x t u u =其中：pdefun 为函数名。
由给定输入变量可计算出,,c f s 这三个函数。
边界条件可以用下面的函数描述(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂+=∂ （4.2） 这样的边值函数可以写为Matlab 函数[,,,](,,,)a a b b x p q p q pdebc x t u u =初始条件数学描述为00(,)u x t u =，编写一个简单的函数即可0()u pdeic x =还可以选择x 和t 的向量，再加上描述这些函数，就可以用pdepe （）函数求解次偏微分方程，需要用如下格式求解(,@,@,@,,)sol pdepe m pdefun pdeic pdebc x t =【例1】 试求下列偏微分方程2111222221220.024()0.17()u u F u u t x u u F u u tx ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩ 其中： 5.7311.46()xx F x e e -=-，且满足初始条件1(,0)1u x =，2(,1)0u x =及边界条件1221(0,)0,(0,)0,(1,)1,(1,)0u u t u t u t t x x∂∂====∂∂解：对照给出的偏微分方程和（4.1），可将原方程改写为111222120.024/()1.*10.17/()u u x F u u u u x F u u t x ∂∂--⎡⎤⎡⎤⎡⎤⎡⎤∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥∂∂-∂∂⎣⎦⎣⎦⎣⎦⎣⎦可知0m =，且1122120.024/()1,,10.17/()u x F u u c f s u x F u u ∂∂--⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥∂∂-⎣⎦⎣⎦⎣⎦编写下面的Matlab 函数function [c,f,s]=c7mpde(x,t,u,du)c=[1;1];y=u(1)-u(2);F=exp(5.73*y)-exp(-11.46*y);s=F*[-1;1]; f=[0.024*du(1);0.17*du(2)];套用（4.2）中的边界条件，可以写出如下的边值方程左边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦，右边界1100.*100u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 编写下面的Matlab 函数function [pa,qa,pb,qb]=c7mpbc(xa,ua,xb,ub,t) pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0,1]; 另外，描述初值的函数function u0=c7mpic(x) u0=[1;0];有了这三个函数，选定x 和t 向量，则可以由下面的程序直接求此微分方程，得出解1u 和2u 。
最新偏微分方程数值解法的MATLAB源码
[ 原创]偏微分方程数值解法的MATLAB源码【更新完毕】说明：由于偏微分的程序都比较长，比其他的算法稍复杂一些，所以另开一贴，专门上传偏微分的程序谢谢大家的支持！其他的数值算法见：..//Announce/Announce.asp?BoardID=209&id=82450041、古典显式格式求解抛物型偏微分方程（一维热传导方程）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 <= 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轴的等分区间数% 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==7C=1;end%计算步长dx=uX/M;%x的步长dt=uT/N;%t的步长x=(0:M)*dx;t=(0:N)*dt;r=C*dt/dx/dx;%步长比r1=1-2*r;if r > 0.5disp('r > 0.5,不稳定')end%计算初值和边值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i));endfor j=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));end%逐层求解for j=1:Nfor i=2:MU(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]=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)=psi2(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]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C);%设置参数C的默认值if nargin==7C=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-2Diag(i)=1+2*r;Low(i)=-r;Up(i)=-r;endDiag(M-1)=1+2*r;%计算初值和边值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i));endfor j=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));end%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）for j=1:Nb1=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,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:线性方程组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:nL(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:nx(i)=b(i)-L(i-1)*x(i-1);end%赶的过程x(n)=x(n)/D(n);for i=n-1:-1:1x(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 <= 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轴的等分区间数%%应用举例：%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-2Diag(i)=1+r;Low(i)=-r/2;Up(i)=-r/2;endDiag(M-1)=1+r;%计算初值和边值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i));endfor j=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));endB=zeros(M-1,M-1);for i=1:M-2B(i,i)=1-r;B(i,i+1)=r/2;B(i+1,i)=r/2;endB(M-1,M-1)=1-r;%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）for j=1:Nb1=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,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) %正方形区域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)% u(x,ub)=psi2(x)%%输出参数：U -解矩阵，第一行表示y=0时的值，第二行表示第y=h时的值……% x -横坐标% y -纵坐标%输入参数：ub -变量边界值的上限% phi1,phi2,psi1,psi2 -边界函数，定义为内联函数% M -横纵坐标的等分区间数% type -求解差分方程的迭代格式，若type='Jacobi'，采用Jacobi迭代格式% 若type='GS'，采用Guass-Seidel迭代格式。
(完整版)偏微分方程的MATLAB解法
引言偏微分方程定解问题有着广泛的应用背景。
人们用偏微分方程来描述、解释或者预见各种自然现象，并用于科学和工程技术的各个领域fll。
然而，对于广大应用工作者来说，从偏微分方程模型出发，使用有限元法或有限差分法求解都要耗费很大的工作量，才能得到数值解。
现在，MATLAB PDEToolbox已实现对于空间二维问题高速、准确的求解过程。
偏微分方程如果一个微分方程中出现的未知函数只含一个自变量，这个方程叫做常微分方程，也简称微分方程；如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。
常用的方法有变分法和有限差分法。
变分法是把定解问题转化成变分问题，再求变分问题的近似解；有限差分法是把定解问题转化成代数方程，然后用计算机进行计算；还有一种更有意义的模拟法，它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。
虽然物理现象本质不同，但是抽象地表示在数学上是同一个定解问题，如研究某个不规则形状的物体里的稳定温度分布问题，由于求解比较困难，可作相应的静电场或稳恒电流场实验研究，测定场中各处的电势，从而也解决了所研究的稳定温度场中的温度分布问题。
随着物理科学所研究的现象在广度和深度两方面的扩展，偏微分方程的应用范围更广泛。
从数学自身的角度看，偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。
从这个角度说，偏微分方程变成了数学的中心。
一、MATLAB方法简介及应用1.1 MATLAB简介MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
1.2 Matlab主要功能数值分析数值和符号计算工程与科学绘图控制系统的设计与仿真数字图像处理数字信号处理通讯系统设计与仿真财务与金融工程1.3 优势特点1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；2) 具有完备的图形处理功能，实现计算结果和编程的可视化；3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。
(完整版)偏微分方程的MATLAB解法
引言偏微分方程定解问题有着广泛的应用背景。
人们用偏微分方程来描述、解释或者预见各种自然现象，并用于科学和工程技术的各个领域fll。
然而，对于广大应用工作者来说，从偏微分方程模型出发，使用有限元法或有限差分法求解都要耗费很大的工作量，才能得到数值解。
现在，MATLAB PDEToolbox已实现对于空间二维问题高速、准确的求解过程。
偏微分方程如果一个微分方程中出现的未知函数只含一个自变量，这个方程叫做常微分方程，也简称微分方程；如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。
常用的方法有变分法和有限差分法。
变分法是把定解问题转化成变分问题，再求变分问题的近似解；有限差分法是把定解问题转化成代数方程，然后用计算机进行计算；还有一种更有意义的模拟法，它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。
虽然物理现象本质不同，但是抽象地表示在数学上是同一个定解问题，如研究某个不规则形状的物体里的稳定温度分布问题，由于求解比较困难，可作相应的静电场或稳恒电流场实验研究，测定场中各处的电势，从而也解决了所研究的稳定温度场中的温度分布问题。
随着物理科学所研究的现象在广度和深度两方面的扩展，偏微分方程的应用范围更广泛。
从数学自身的角度看，偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。
从这个角度说，偏微分方程变成了数学的中心。
一、MATLAB方法简介及应用1.1 MATLAB简介MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
1.2 Matlab主要功能数值分析数值和符号计算工程与科学绘图控制系统的设计与仿真数字图像处理数字信号处理通讯系统设计与仿真财务与金融工程1.3 优势特点1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；2) 具有完备的图形处理功能，实现计算结果和编程的可视化；3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。
Matlab求解微分方程(组)及偏微分方程(组)
第四道 Matlab供解微分圆程（组）之阳早格格创做表里介绍：Matlab供解微分圆程（组）下令供解真例：Matlab供解微分圆程（组）真例本质应用问题通过数教修模所归纳得到的圆程，绝大普遍皆是微分圆程，真真能得到代数圆程的机会很少.另一圆里，不妨供解的微分圆程也是格中有限的，特地是下阶圆程战偏偏微分圆程（组）.那便央供咱们必须钻研微分圆程（组）的解法：剖析解法战数值解法.一．相关函数、下令及简介1.正在Matlab中，用大写字母D表示导数，Dy表示y关于自变量的一阶导数，D2y表示y关于自变量的二阶导数，依此类推.函数dsolve用去办理常微分圆程（组）的供解问题，调用要领为：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve用去解标记常微分圆程、圆程组，如果不初初条件，则供出通解，如果有初初条件，则供出特解.注意，系统缺省的自变量为t2.函数dsolve供解的是常微分圆程的透彻解法，也称为常微分圆程的标记解.然而是，有洪量的常微分圆程虽然从表里上道，其解是存留的，然而咱们却无法供出其剖析解，此时，咱们需要觅供圆程的数值解，正在供常微分圆程数值解圆里，MATLAB 具备歉富的函数，咱们将其统称为solver ，其普遍要领为：[T,Y]=solver(odefun,tspan,y0)证明：(1)solver 为下令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一.(2)odefun 是隐现微分圆程'(,)y f t y =正在积分区间tspan 0[,]f t t =上从0t 到f t 用初初条件0y 供解.(3)如果要赢得微分圆程问题正在其余指定时间面012,,,,f t t t t 上的解，则令tspan 012[,,,]f t t t t =(央供是单调的).(4)果为不一种算法不妨灵验的办理所有的ODE 问题，为此，Matlab 提供了多种供解器solver ，对付于分歧的ODE 问题，采与分歧的solver.表1 Matlab 华文本文献读写函数证明：ode23、ode45是极其时常使用的用去供解非刚刚性的尺度形式的一阶微分圆程（组）的初值问题的解的Matlab时常使用步调，其中：ode23采与龙格-库塔2阶算法，用3阶公式做缺面预计去安排步少，具备矮等的粗度.ode45则采与龙格-库塔4阶算法，用5阶公式做缺面预计去安排步少，具备中等的粗度.3．正在matlab下令窗心、步调或者函数中创修局部函数时，可用内联函数inline，inline函数形式相称于编写M 函数文献，然而不需编写M-文献便不妨形貌出某种数教关系.调用inline函数，只可由一个matlab表白式组成，而且只可返回一个变量，不允许[u,v]那种背量形式.果而，所有央供逻辑运算或者乘法运算以供得最后截止的场合，皆不克不迭应用inline函数，inline函数的普遍形式为：FunctionName=inline(‘函数真质’, ‘所有自变量列表’)比圆：（供解F(x)=x^2*cos(a*x)-b ,a,b是标量；x是背量）正在下令窗心输进:Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’);g= Fofx([pi/3 pi/3.5],4,1)注意：由于使用内联对付象函数inline 不需要其余修坐m 文献，所有使用比较便当，其余正在使用ode45函数的时间，定义函数往往需要编写一个m 文献去单独定义，那样便当于管造文献，那里不妨使用inline 去定义函数.二．真例介绍例1 供解微分圆程2'2x y xy xe -+= 步调：syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x ’) 例2 供微分圆程'0x xy y e +-=正在初初条件(1)2y e =下的特解并画出解函数的图形.步调：syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x ’);ezplot(y)例 3 供解微分圆程组530t dx x y e dt dy x y dt ⎧++=⎪⎪⎨⎪--=⎪⎩正在初初条件00|1,|0t t x y ====下的特解并画出解函数的图形.步调：syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')simple(x);simple(y)ezplot(x,y,[0,1.3]);axis auto2.用ode23、ode45等供解非刚刚性尺度形式的一阶微分圆程（组）的初值问题的数值解（近似解）例4 供解微分圆程初值问题2222(0)1dy y x x dx y ⎧=-++⎪⎨⎪=⎩的数值解，供解范畴为区间[0,0.5].步调：fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1);plot(x,y,'o-')例5 供解微分圆程22'2(1)0,(0)1,(0)0d y dy y y y y dt dt μ--+===的解，并画出解的图形.12,,7dy x y x dt μ===，则编写M-文献function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)];end正在Matlab 下令窗心编写步调y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或者[t,x]=ode45('vdp',[0,40],y0);y=x(:,1);dy=x(:,2);plot(t,y,t,dy)训练与思索：M-文献vdp.m 改写成inline 函数步调？Euler 合线法供解的基础思维是将微分圆程初值问题 化成一个代数(好分)圆程，主要步调是用好商()()y x h y x h +-代替微商dydx ，于是 记1,(),k k k k x x h y y x +=+=进而1(),k k y y x h +=+于是例6用Euler 合线法供解微分圆程初值问题的数值解（步少h 与），供解范畴为区间[0,2].分解：本问题的好分圆程为步调：>> clear>> f=sym('y+2*x/y^2');>> a=0;>> b=2;>> h=0.4;>> n=(b-a)/h+1;>> x=0;>> y=1;>> szj=[x,y];%数值解>> for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});%subs ，替换函数x=x+h;szj=[szj;x,y];end>>szj>> plot(szj(:,1),szj(:,2))证明：替换函数subs比圆：输进subs(a+b,a,4) 意义便是把a用4替换掉，返回4+b，也不妨替换多个变量，比圆：subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分别用字符alpha替换a战2替换b，返回 cos(alpha)+sin(2)特地证明：本问题可进一步利用四阶Runge-Kutta法供解，Euler合线法本质上便是一阶Runge-Kutta法，Runge-Kutta法的迭代公式为相映的Matlab步调为：>> clear>> f=sym('y+2*x/y^2');>> a=0;>> b=2;>> h=0.4;>> n=(b-a)/h+1;>> x=0;>> y=1;>> szj=[x,y];%数值解>> for i=1:n-1l1=subs(f,{'x','y'},{x,y});替换函数l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2});l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2});l4=subs(f,{'x','y'},{x+h,y+l3*h});y=y+h*(l1+2*l2+2*l3+l4)/6;x=x+h;szj=[szj;x,y];end>>szj>> plot(szj(:,1),szj(:,2))训练与思索：(1)ode45供解问题并比较好别.(2)利用Matlab 供微分圆程(4)(3)''20y y y -+=的解. (3)供解微分圆程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解.(4)利用Matlab 供微分圆程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解.指示：尽大概多的思量解法三．微分圆程变更为一阶隐式微分圆程组Matlab 微分圆程解算器只可供解尺度形式的一阶隐式微分圆程（组）问题，果此正在使用ODE 解算器之前，咱们需要干的第一步，也是最要害的一步便是借帮状态变量将微分圆程（组）化成Matlab 可交受的尺度形式.天然，如果ODEs 由一个或者多个下阶微分圆程给出，则咱们应先将它变更成一阶隐式常微分圆程组.底下咱们以二个下阶微分圆程组形成的ODEs 为例介绍怎么样将它变更成一个一阶隐式微分圆程组.Step 1 将微分圆程的最下阶变量移到等式左边，其余移到左边，并按阶次从矮到下排列.形式为：Step 2 为每一阶微分式采用状态变量，最下阶除中注意：ODEs 中所有是果变量的最下阶次之战便是需要的状态变量的个数，最下阶的微分式不需要给它状态变量.Step 3 根据采用的状态变量，写出所有状态变量的一阶微分表白式训练与思索：(1)供解微分圆程组 其中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(2)供解隐式微分圆程组提示：使用标记预计函数solve供'''',x y ，而后利用供解微分圆程的要领四．偏偏微分圆程解法Matlab 提供了二种要领办理PDE 问题，一是使用pdepe 函数，它不妨供解普遍的PDEs,具备较大的通用性，然而只收援下令形式调用；二是使用PDE 工具箱，不妨供解特殊PDE 问题，PDEtoll 有较大的限造性，比圆只可供解二阶PDE 问题，而且不克不迭办理片微分圆程组，然而是它提供了GUI 界里，从搀纯的编程中解脱出去，共时还不妨通过File —>Save As 直交死成M 代码.1.普遍偏偏微分圆程（组）的供解(1)Matlab 提供的pdepe 函数，不妨直交供解普遍偏偏微分圆程（组），它的调用要领为：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题形貌函数，它必须换成尺度形式：那样，PDE 便不妨编写出心函数：[c,f,s]=pdefun(x,t,u,du)，m,x,t 对付应于式中相关参数，du 是u 的一阶导数，由给定的输进变量可表示出c,f,s 那三个函数.@pdebc 是PDE 的鸿沟条件形貌函数，它必须化为形式：于是边值条件不妨编写函数形貌为：[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a 表示下鸿沟，b 表示上鸿沟.@pdeic 是PDE 的初值条件，必须化为形式：00(,)u x t u =，故不妨使用函数形貌为：u0=pdeic(x)sol 是一个三维数组，sol(:,:,i)表示i u 的解，换句话道，k u 对付应x(i)战t(j)时的解为sol(i,j,k)，通过sol ，咱们不妨使用pdeval 函数直交预计某个面的函数值.(2)真例证明供解偏偏微分其中， 5.7311.46()x x F x e e -=-且谦脚初初条件12(,0)1,(,0)0u x u x ==及鸿沟条件1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0u u t u t t x ∂===∂解：(1)对付照给出的偏偏微分圆程战pdepe 函数供解的尺度形式，本圆程改写为 可睹1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦%目标PDE 函数function [c,f,s]=pdefun(x,t,u,du)c=[1;1];f=[0.024*du(1);0.17*du(2)];temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp))end(2)鸿沟条件改写为：下鸿沟2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上鸿沟1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ %鸿沟条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t)pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0;1];end(3)初值条件改写为：1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值条件函数function u0=pdeic(x)u0=[1;0];end(4)编写主调函数clcx=0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);subplot(2,1,1)surf(x,t,sol(:,:,1))subplot(2,1,2)surf(x,t,sol(:,:,2))训练与思索： This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(1)PDEtool （GUI ）供解偏偏微分圆程的普遍步调正在Matlab 下令窗心输进pdetool ，回车，PDE 工具箱的图形用户界里(GUI)系统便开用了.从定义一个偏偏微分圆程问题到完毕解偏偏微分圆程的定解，所有历程大概不妨分为六个阶段Step 1 “Draw模式”画造仄里有界地区Ω，通过公式把Matlab系统提供的真体模型：矩形、圆、椭圆战多边形，拉拢起去，死成需要的仄里地区.Step 2 “Boundary模式”定义鸿沟，声明分歧鸿沟段的鸿沟条件.Step 3 “PDE模式”定义偏偏微分圆程，决定圆程典型战圆程系数c,a,f,d，根据简直情况，还不妨正在分歧子地区声明分歧系数.Step 4 “Mesh模式”网格化地区Ω，不妨统造自动死成网格的参数，对付死成的网格举止多次细化，使网格分隔更细更合理.Step 5 “Solve模式”解偏偏微分圆程，对付于椭圆型圆程不妨激活并统造非线性自符合解题器去处理非线性圆程；对付于扔物线型圆程战单直型圆程，树坐初初鸿沟条件后不妨供出给定时刻t的解；对付于特性值问题，不妨供出给定区间上的特性值.供解完毕后，不妨返回到Step 4，对付网格进一步细化，举止再次供解.Step 6 “View模式”预计截止的可视化，不妨通过树坐系统提供的对付话框，隐现所供的解的表面图、网格图、等下线图战箭头梯形图.对付于扔物线型战单直线型问题的解还不妨举止径画演示.(2)真例证明用法供解一个正圆形地区上的特性值问题：正圆形地区为：11,1 1.-≤≤-≤≤x x(1)使用PDE工具箱挨开GUI供解圆程(2)加进Draw模式，画造一个矩形，而后单打矩形，正在弹出的对付话框中树坐Left=-1,Bottom=-1,Width=2,Height=2,确认并关关对付话框(3)加进Boundary模式，鸿沟条件采与Dirichlet条件的默认值(4)加进PDE模式，单打工具栏PDE按钮，正在弹出的对付话框中圆程典型采用Eigenmodes,参数树坐c=1,a=-1/2,d=1,确认后关关对付话框(5)单打工具栏的∆按钮，对付正圆形地区举止初初网格剖分，而后再对付网格进一步细化剖分一次(6)面开solve菜单，单打Parameters选项，正在弹出的对付话框中树坐特性值地区为[-20,20](7)单打Plot菜单的Parameters项，正在弹出的对付话框中选中Color、Height(3-D plot)战show mesh项，而后单打Done确认(8)单打工具栏的“=”按钮，开初供解。
偏微分的MATLAB求解
MATLAB提供了dpepe函数来求解该问题的数值解。其基本调用格式为： sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
2.1 用偏微分方程工具箱求解微分方程
（5）求解此双曲问题 u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d); %得到如下结果： %Time: 0.166667 %Time: 0.333333 %Time: 0.5 %… %Time: 4.5 %Time: 4.66667 %Time: 4.83333 %Time: 5 %428 successful steps %62 failed attempts %982 function evaluations %1 partial derivatives %142 LU decompositions %981 solutions of linear systems %现在把解u1用动态图形表示出来． （6）矩形网格插值 delta=-1:0.1:1; [uxy, tn, a2, a3]=tri2grid(p, t, u1(:, 1), delta, delta); gp=[tn; a2; a3];
29图2210a图2210b仿真解与精确解的误差图30223定解问题的计算机仿真本小节主要对定解问题的求解给出计算机仿真显示直观明了地给出解动态或静态分布有利于对物理模型和物理意义的理解
偏微分与 MATLAB
本章将主要讲述如何用MATLAB实现对偏微分 方程的仿真求解 ．MATLAB的偏微分方程工具箱 （PDE Toolbox）的出现，为偏微分方程的求解以 及定性研究提供了捷径．主要步骤为：
【精品】偏微分的MATLAB数值解法课件
方法一：pdepe函数实现
• x=0:1:40; • t=0:0.01:0.2; • m=0; • sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t); • b=sol(20,:); • plot(x,b); • title('the solution of u') • xlabel('x') • ylabel('y') • zlabel('u')
偏微分的MATLAB数值解法
偏微分的MATLAB数值解法
• 方法一：pdepe函数实现 • 方法二：pdetool实现 • 方法三：程序实现
方法一：pdepe函数实现
• @pdeic： • function u0=pdeic(x) • if x<10 • u0=0; • elseif x<30 • u0=1; • else • u0=0; • end
•
end
• end
方法三：程序实现
图 22.12 波动方程解析解的分布
偏微分的MATLAB数值解法
• 方法总结： • 1.pdede调用简单，但计算功能稍弱 • 2.pdetool使用方便，但限于四种方程类
型 • 3.程序编写较为繁琐
方法一：pdepe函数实现
方法二：pdetool实现
• 1.pdetool界面 • 2.选定求解微分方程类型（双曲线、抛物线、椭
圆、特殊值型）并设定参数 • 3.绘制求解区域 • 4.边界条件和初值条件（Dirichlet和Neumann） • 5.生成网格 • 6.求解方程并绘制图形
方法二：pdetool实现
• 应用实例：
u(ux,y)
x2 y x0
偏微分方程(PDEs)的MATLAB数值解法
偏微分方程的MATLAB求解精讲©MA TLAB求解微分/偏微分方程，一直是一个头大的问题，两个字，“难过”，由于MA TLAB对LaTeX的支持有限，所有方程必须化成MA TLAB可接受的标准形式，不支持像其他三个数学软件那样直接傻瓜式输入，这个真把人给累坏了！不抱怨了，还是言归正传，回归我们今天的主体吧！MA TLAB提供了两种方法解决PDE问题，一是pdepe()函数，它可以求解一般的PDEs，据用较大的通用性，但只支持命令行形式调用。
二是PDE工具箱，可以求解特殊PDE问题，PDEtool有较大的局限性，比如只能求解二阶PDE问题，并且不能解决偏微分方程组，但是它提供了GUI界面，从繁杂的编程中解脱出来了，同时还可以通过File->Save As直接生成M代码一、一般偏微分方程组(PDEs)的MA TLAB求解 (3)1、pdepe函数说明 (3)2、实例讲解 (4)二、PDEtool求解特殊PDE问题 (6)1、典型偏微分方程的描述 (6)（1）椭圆型 (6)（2）抛物线型 (6)（3）双曲线型 (6)（4）特征值型 (7)2、偏微分方程边界条件的描述 (8)（1）Dirichlet条件 (8)（2）Neumann条件 (8)3、求解实例 (9)一、一般偏微分方程组(PDEs)的MATLAB 求解1、pdepe 函数说明MA TLAB 语言提供了pdepe()函数，可以直接求解一般偏微分方程(组)，它的调用格式为sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)【输入参数】@pdefun ：是PDE 的问题描述函数，它必须换成下面的标准形式(,,)[(,,,)](,,,)()m m u u u uc x t x x f x t u s x t u x t x x x−∂∂∂∂∂=+∂∂∂∂∂式1 这样，PDE 就可以编写下面的入口函数 [c,f,s]=pdefun(x,t,u,du)m,x,t 就是对应于(式1)中相关参数，du 是u 的一阶导数，由给定的输入变量即可表示出出c,f,s 这三个函数@pdebc ：是PDE 的边界条件描述函数，必须先化为下面的形式(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂+=∂ 于是边值条件可以编写下面函数描述为 [pa,qa,pb,qb]=pdebc(x,t,u,du)其中a 表示下边界，b 表示下边界@pdeic ：是PDE 的初值条件，必须化为下面的形式00(,)u x t u =我们使用下面的简单的函数来描述为 u0=pdeic(x)m,x,t ：就是对应于(式1)中相关参数【输出参数】sol ：是一个三维数组，sol(:,:,i)表示u i 的解，换句话说u k 对应x(i)和t(j)时的解为sol(i,j,k)通过sol ，我们可以使用pdeval()直接计算某个点的函数值2、实例讲解试求解下面的偏微分2111222221220.024()0.17()u u F u u t xu u F u u tx ∂∂=−− ∂∂ ∂∂ =−− ∂∂ 其中， 5.7311.46()x x F x e e −=−，且满足初始条件12(,0)1,(,0)0u x u x ==及边界条件1221(0,)0,(0,)0,(1,)1,(1,)0u ut u t u t t x x∂∂====∂∂【解】（1）对照给出的偏微分方程，根据标注形式，则原方程可以改写为111222120.024()1.*1()0.17u u F u u x u u F u u t t x ∂−−∂∂∂=+ ∂−∂∂∂可见m=0，且1122120.024()1,,1()0.17u F u u x c f s u F u u x ∂−− ∂===∂−∂%% 目标PDE 函数function [c,f,s]=pdefun (x,t,u,du) c=[1;1];f=[0.024*du(1);0.17*du(2)]; temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp));（2）边界条件改写为12011010.*.*00000u f f u −+=+=下边界上边界%% 边界条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t) %a 表示下边界，b 表示上边界 pa=[0;ua(2)];qa=[1;0]; pb=[ub(1)-1;0]; qb=[0;1];（3）初值条件改写为1210u u =%% 初值条件函数function u0=pdeic(x) u0=[1;0];（4）最后编写主调函数 clcx=0:0.05:1; t=0:0.05:2; m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);figure('numbertitle','off','name','PDE Demo ——by Matlabsky') subplot(211)surf(x,t,sol(:,:,1)) title('The Solution of u_1') xlabel('X') ylabel('T') zlabel('U') subplot(212)surf(x,t,sol(:,:,2)) title('The Solution of u_2') xlabel('X') ylabel('T') zlabel('U')二、PDEtool 求解特殊PDE 问题MATLAB 的偏微分工具箱(PDE toolbox)可以比较规范的求解各种常见的二阶偏微分方程，但是惋惜的是只能求解特殊二阶的PDE 问题，并且不支持偏微分方程组！PDE toolbox 支持命令行形式求解PDE 问题，但是要记住那些命令以及调用形式真的很累人，还好MATLAB 提供了GUI 可视交互界面pdetool ，在pdetool 中可以很方便的求解一个PDE 问题，并且可以帮我们直接生成M 代码(File->Save As)。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
求解双曲型方程的例子
例24.2.1 用 MATLAB 求解下面波动方程定解问题并动态显示解的分布
2u (2u t 2 x2
2u ) 0 y 2
u
|x
1
u
|x1
0,
u y
y 1
u y
y1 0
π
π
u(x,
y, 0)
atan[ sin(
2
x)], ut
( x,
y,
0)
2
cos(πx)
保持在100 °C，板的右边热量从板向环境空气定常流动，
t t 其他边及内孔边界保持绝缘。初始
°C ，于是概括为如下定解问题;
是板的温度为0 0
d u u0 , t
u 1 0 0 ,在 左 边 界 上
u 1， 在 右 边 界 上 n u = 0， 其 他 边 界 上 n
u t to 0
区域的边界顶点坐标为(-0.5,-0.8), (0.5,-0.8), (-0.5,0.8), (0.5,0.8)。 内边界顶点坐标(-0.05,-0.4), (-0.05,0.4) ,(0.05,-0.4), (0.05,0.4)。
第七步：单击Plot菜单中Parameter选项，打开Plot Selection对话框，选中Color,Height(3D plot)和 Show mesh三项.再单击Polt按钮，显示三维图形解， 如图22.5所示.
第八步：若要画等值线图和矢量场图，单击plot菜单 中parameter 选项，在plot selection对话框中选中 contour 和arrow两选项。然后单击plot按钮，可显示 解的等值线图和矢量场图，如图2.6所示。
网格划分，细化
Solve，Plot
如果有初始条件（与t有关），则在 Solve的Parameters里有其设定，如果 没有初始条件（与t无关），则不必设 定Plot只是确定画图的参数，包括是否 动画，是否3D，是否画出等温线，是否 有箭头。。。
Save As
保存成M-file，自动生成
例： 解热传导方程 ut u f
第五步：选择Mesh菜单中Initialize Mesh命使网格密集化，
如图22.3．
图 22.3 网格密集化
第六步： 解偏微分方程并显示图形解
选择Solve菜单中Solve PDE命令，解 偏微分方程并显示图形解，如图 2.4 所示
例如，对于细杆导热，虽然是一维问题， 可以将宽度y虚拟出来，对应于y的边界 条件和初始条件按照题意制定
Boundary Mode
PDE Mode
PDE Specification，确定偏 微分方程类型共有四种：
椭圆形Elliptic
抛物型Parabolic
双曲型Hyperbolic
Mesh Mode
图 2.6 解的等值线图和矢量场图
求解椭圆型方程的例子
单 位 圆 上 的 poisson方 程 边 值 问 题 ：
-u =1 , = ( x, y ) x 2 y 2 1 ,
u 0 问题的精确解为
u(x,y)= (1 x 2 y 2 ) . 4
求解抛物型方程的例子
考虑一个带有矩形孔的金属板上的热传导问题。板的左边
第三步：选取边界
首先选择Boundary菜单中Boundary Mode 命令，进入边界模式．然后单击Boundary菜单 中Remove All Subdomain Borders选项。从而 去掉子域边界，如图22 2．单击Boundary菜单 中Specify Boundary Conditions选项，打开 Boundary Conditions对话框，输入边界件．本 例取默认条件，即将全部边界设为齐次Dirichlet 条件，边界显示为红色． 如果想将几何与边界信息存储，可选Boundary 菜单中的Export Decomposed Geometrv．Boundary Cond's命令，将它们分 别存储在g、b变量中，并通过MATLAB形成M文 件．
偏微分方程的matlab解法
PDEToolbox注意事项
只能解决二维模型，一维的扩成二维，三 维的缩成二维，时间维不计算在内 公式类型，只能解决部分偏微分方程，由 公式类型决定 边界条件两种，Dirichlet和Neumann 初始条件
先确定方程大类
Draw Mode
画图模式，先将处理的区域画出来，二 维，方形，圆形，支持多边形，可以手 动更改坐标，旋转rotate
图 22.2 定解问题的边界
第四步：设置方程类型
选择PDE菜单中PDE Mode命令，进入PDE模式， 再单击PDE菜单中PDE Secification选项，打开 PDE Secification对话框，设置方程类型．
本例取抛物型方程 du(cu)auf,
t
故参数c,a,f,d,分别是l，0，10，1．
exp[cos(
2
y)]
已知求解域是方形区域，空间坐标的个数由具体问题确定．
此课件下载可自行编辑修改，仅供参考！ 感谢您的支持，我们努力做得更好！谢谢
边界条件是齐次类型，定解区域自定。
【解】 第一步：启动MATLAB，键入命令pdetool并回车，就进 入GUI．在Options菜单下选择Gid命令，打开栅格，栅 格使用户容易确定所绘图形的大小． 第二步：选定定解区域本题为自定区域：自拟定解区 域如图22 1所示：E1-E2+R1-E3．具体用快捷工具分别 画椭圆E1、圆E2、矩形R1、圆E3．然后在Set formula 栏中进行编辑并用算术运算符将图形对象名称连接起 来(或删去默认的表达式，直接键入E1-E2+R1-E3)