一维偏微分方程的pdepe(matlab)函数解法

本文根据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.


相关文章

  • 线性谐振子的不同解法比较
  • 线性谐振子的不同解法比较 关键词:一维谐振子:能量本征值:波函数 摘 要:一维线性谐振子作为量子力学中的基础模型,它的解决方法具有多样性并随着科学工作者的努力和对数学理论的应用的不断深入(如群论和群表示理论),谐振子的解法将会最优化,并会对 ...查看


  • Matlab中常微分方程数值解法
  • Matlab中常微分方程数值解法讲解 作者:dynamic 时间:2008.12.10 版权:All Rights Reserved By www.matlabsky.cn ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ ...查看


  • 算法大全第15章常微分方程的解法
  • 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步,通常需要求出方程的解来说明实际现象,并加以检验.如果能得到解析形式的解固然是便于分析和应用的,但是我们知道,只有线性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以肯定得 ...查看


  • 矩阵分解及无约束最优化方法
  • 矩阵分解及无约束最优化 方法的原理和应用简介 --最优化方法课程实验报告 学 院:数学与统计学院 班 级:硕2041班 姓 名:王彭 学 号:指导教师:阮小娥 同 组 人:陈莹 钱东东 矩阵分解及无约束最优化方法的原理和应用简介 矩阵分解及 ...查看


  • 数值分析学习心得体会
  • 数值分析学习感想 一个学期的数值分析,在老师的带领下,让我对这门课程有了深刻的理解和感悟.这门 课程是一个十分重视算法和原理的学科,同时它能够将人的思维引入数学思考的模式,在处 理问题的时候,可以合理适当的提出方案和假设.他的内容贴近实际, ...查看


  • matlab中有关数学建模的函数
  • intersect(a,b):返回向量a ,b 的公共部分 intersect (A,B,'row' ):A,B为相同列数的矩阵,返回元素相同的行 [c,ia,ib]=intersect(a,b):c为a ,b 的公共元素,ia 表示公共元 ...查看


  • 常微分方程常用数值解法
  • 第一章 绪论 1.1 引言 常微分方程是现代数学的一个重要分支,是人们解决各种实际问题的有效工具.微分方程的理论和方法从17世纪末开始发展起来,很快成了研究自然现象的强有力工具,在17到18世纪,在力学.天文.科学技术.物理中,就已借助微分 ...查看


  • 数值分析作业思考题
  • 数值分析思考题1 1.讨论绝对误差(限).相对误差(限)与有效数字之间的关系. 2.相对误差在什么情况下可以用下式代替? e *x *-x e =*= x x * *r 3.查阅何谓问题的"病态性",并区分与" ...查看


  • 二维热传导方程有限容积法的MATLAB实现_薛琼
  • ComputerEngineeringandApplications计算机工程与应用2012,48(24)197 二维热传导方程有限容积法的MATLAB实现 薛琼1,肖小峰2 XUEQiong1,XIAOXiaofeng2 1.武汉理工大学 ...查看


热门内容