本文根据matlab 帮助进行加工,根据
matlab 帮助上的例子,帮助更好的理解一维偏微分方程的pdepe 函数解法,主要加工在于程序的注释上。
Examples
Example 1.
This example illustrates the straightforward
formulation, computation, and plotting of the solution of a single PDE.
This equation holds on an interval The PDE satisfies the initial condition
and boundary conditions
for times
.
It is convenient to use subfunctions to place all the functions required by pdepe in a single function. function pdex1
m = 0;
x = linspace(0,1,20);
%linspace(x1,x2,N)linspace是Matlab 中的一个指令,用于产生x1,x2之间的N 点行矢量。
%其中x1、x2、N 分别为起始值、终止值、元素个数。若缺省N ,默认点数为100 t = linspace(0,2,5);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % Extract the first solution component as u. u = sol(:,:,1);
% A surface plot is often a good way to study a solution. surf(x,t,u)
title('Numerical solution computed with 20 mesh points.') xlabel('Distance x') ylabel('Time t')
% A solution profile can also be illuminating. figure
plot(x,u(end,:))
title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)')
% -------------------------------------------------------------- function [c,f,s] = pdex1pde(x,t,u,DuDx) c = pi^2; f = DuDx; s = 0;
% -------------------------------------------------------------- function u0 = pdex1ic(x) u0 = sin(pi*x);
% -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0;
pr = pi * exp(-t); qr = 1;
In this example, the PDE, initial condition, and boundary conditions are coded in subfunctions pdex1pde , pdex1ic , and pdex1bc . The surface plot shows the behavior of the solution.
The following plot shows the solution profile at the final value of t (i.e., t = 2) .
2
f (u ) =u 我们再将该问题复杂化,比如在原方程右边加一项,
对于标准形式
s =u 2,其余条件不变
function pdex1
m = 0;
x = linspace(0,1,20);
%linspace(x1,x2,N)linspace是Matlab 中的一个指令,用于产生x1,x2之间的N 点行矢量。
%其中x1、x2、N 分别为起始值、终止值、元素个数。若缺省N ,默认点数为100 t = linspace(0,2,5);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % Extract the first solution component as u. u = sol(:,:,1);
% A surface plot is often a good way to study a solution. surf(x,t,u)
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x') ylabel('Time t')
% A solution profile can also be illuminating. figure
plot(x,u(end,:))
title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)')
% -------------------------------------------------------------- function [c,f,s] = pdex1pde(x,t,u,DuDx) c = pi^2; f = DuDx; s = u^2;
% -------------------------------------------------------------- function u0 = pdex1ic(x) u0 = sin(pi*x);
% -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0;
pr = pi * exp(-t); qr = 1;
对比结果有差异,表示可行
Example 2. This
example illustrates the solution of a system of
PDEs. The
problem has boundary layers at both ends of the interval. The solution changes rapidly for small
. The PDEs are
where
This equation holds on an interval The PDE satisfies the initial conditions
.
for times
.
and boundary conditions
In the form expected by pdepe
, the equations are
The boundary conditions on the partial derivatives of have to be written in terms of the flux. In the form expected by
pdepe , the left boundary condition is
and the right boundary condition is
The solution changes rapidly for small . The program selects the step size in time to resolve this sharp change, but to see this behavior in the plots, the example must select the output times accordingly. There are boundary layers in the solution at both ends of [0,1], so the example places mesh points near 0 and 1 to resolve these sharp changes. Often some experimentation is needed to select a mesh that reveals the behavior of the solution. 程序在pdex4.m 中
function pdex4 %pdex4为文件名 clc m = 0;
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];
sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
%pdepe函数,用于直接求解偏微分方程,其形式为sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
%下面为作图步骤 u1 = sol(:,:,1);
% ui = sol(:,:,i) is an approximation to the ith component of the
solution vector u . u2 = sol(:,:,2);
figure
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x') ylabel('Time t')
figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x') ylabel('Time t')
% -------------------------------------------------------------- function [c,f,s] = pdex4pde(x,t,u,DuDx)
%被调用的函数,其作用是描述偏微分方程 c = [1; 1];
f = [0.024; 0.17] .* DuDx; y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y); s = [-F; F];
% -------------------------------------------------------------- function u0 = pdex4ic(x);
%此函数是用来描述初值 u0 = [1; 0];
% -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t) %此函数用来描述边界条件 pl = [0; ul(2)]; ql = [1; 0];
pr = [ur(1)-1; 0]; qr = [0; 1];
In this example, the PDEs, initial conditions, and boundary conditions are coded in subfunctions pdex4pde , pdex4ic , and pdex4bc .
The surface plots show the behavior of the solution components.
本文根据matlab 帮助进行加工,根据
matlab 帮助上的例子,帮助更好的理解一维偏微分方程的pdepe 函数解法,主要加工在于程序的注释上。
Examples
Example 1.
This example illustrates the straightforward
formulation, computation, and plotting of the solution of a single PDE.
This equation holds on an interval The PDE satisfies the initial condition
and boundary conditions
for times
.
It is convenient to use subfunctions to place all the functions required by pdepe in a single function. function pdex1
m = 0;
x = linspace(0,1,20);
%linspace(x1,x2,N)linspace是Matlab 中的一个指令,用于产生x1,x2之间的N 点行矢量。
%其中x1、x2、N 分别为起始值、终止值、元素个数。若缺省N ,默认点数为100 t = linspace(0,2,5);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % Extract the first solution component as u. u = sol(:,:,1);
% A surface plot is often a good way to study a solution. surf(x,t,u)
title('Numerical solution computed with 20 mesh points.') xlabel('Distance x') ylabel('Time t')
% A solution profile can also be illuminating. figure
plot(x,u(end,:))
title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)')
% -------------------------------------------------------------- function [c,f,s] = pdex1pde(x,t,u,DuDx) c = pi^2; f = DuDx; s = 0;
% -------------------------------------------------------------- function u0 = pdex1ic(x) u0 = sin(pi*x);
% -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0;
pr = pi * exp(-t); qr = 1;
In this example, the PDE, initial condition, and boundary conditions are coded in subfunctions pdex1pde , pdex1ic , and pdex1bc . The surface plot shows the behavior of the solution.
The following plot shows the solution profile at the final value of t (i.e., t = 2) .
2
f (u ) =u 我们再将该问题复杂化,比如在原方程右边加一项,
对于标准形式
s =u 2,其余条件不变
function pdex1
m = 0;
x = linspace(0,1,20);
%linspace(x1,x2,N)linspace是Matlab 中的一个指令,用于产生x1,x2之间的N 点行矢量。
%其中x1、x2、N 分别为起始值、终止值、元素个数。若缺省N ,默认点数为100 t = linspace(0,2,5);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t); % Extract the first solution component as u. u = sol(:,:,1);
% A surface plot is often a good way to study a solution. surf(x,t,u)
title('Numerical solution computed with 20 mesh points.')
xlabel('Distance x') ylabel('Time t')
% A solution profile can also be illuminating. figure
plot(x,u(end,:))
title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)')
% -------------------------------------------------------------- function [c,f,s] = pdex1pde(x,t,u,DuDx) c = pi^2; f = DuDx; s = u^2;
% -------------------------------------------------------------- function u0 = pdex1ic(x) u0 = sin(pi*x);
% -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0;
pr = pi * exp(-t); qr = 1;
对比结果有差异,表示可行
Example 2. This
example illustrates the solution of a system of
PDEs. The
problem has boundary layers at both ends of the interval. The solution changes rapidly for small
. The PDEs are
where
This equation holds on an interval The PDE satisfies the initial conditions
.
for times
.
and boundary conditions
In the form expected by pdepe
, the equations are
The boundary conditions on the partial derivatives of have to be written in terms of the flux. In the form expected by
pdepe , the left boundary condition is
and the right boundary condition is
The solution changes rapidly for small . The program selects the step size in time to resolve this sharp change, but to see this behavior in the plots, the example must select the output times accordingly. There are boundary layers in the solution at both ends of [0,1], so the example places mesh points near 0 and 1 to resolve these sharp changes. Often some experimentation is needed to select a mesh that reveals the behavior of the solution. 程序在pdex4.m 中
function pdex4 %pdex4为文件名 clc m = 0;
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];
sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
%pdepe函数,用于直接求解偏微分方程,其形式为sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
%下面为作图步骤 u1 = sol(:,:,1);
% ui = sol(:,:,i) is an approximation to the ith component of the
solution vector u . u2 = sol(:,:,2);
figure
surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x') ylabel('Time t')
figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x') ylabel('Time t')
% -------------------------------------------------------------- function [c,f,s] = pdex4pde(x,t,u,DuDx)
%被调用的函数,其作用是描述偏微分方程 c = [1; 1];
f = [0.024; 0.17] .* DuDx; y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y); s = [-F; F];
% -------------------------------------------------------------- function u0 = pdex4ic(x);
%此函数是用来描述初值 u0 = [1; 0];
% -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t) %此函数用来描述边界条件 pl = [0; ul(2)]; ql = [1; 0];
pr = [ur(1)-1; 0]; qr = [0; 1];
In this example, the PDEs, initial conditions, and boundary conditions are coded in subfunctions pdex4pde , pdex4ic , and pdex4bc .
The surface plots show the behavior of the solution components.