一维抛物线偏微分方程数值 解法(3)(附图及matlab程序)
一维抛物型偏微分方程初边值问题求解
一维抛物型偏微分方程初边值问题求解摘要:一、引言二、一维抛物型偏微分方程1.定义与性质2.初边值问题三、求解方法1.紧差分格式2.追赶法3.有限元算法四、Matlab程序实现1.紧差分格式程序2.追赶法程序五、结论与展望正文:一、引言在数学、物理等领域,偏微分方程是一类重要的方程。
其中,一维抛物型偏微分方程在科学研究和实际应用中具有广泛的意义。
本文将探讨一维抛物型偏微分方程的初边值问题的求解方法,并介绍相应的Matlab程序实现。
二、一维抛物型偏微分方程1.定义与性质一维抛物型偏微分方程是指具有如下形式的方程:u_t = a * u_xx其中,u(x, t) 表示未知函数,t 表示时间,x 表示空间坐标,a 为常数。
2.初边值问题初边值问题是指在给定的初始条件和边界条件下求解偏微分方程的问题。
在一维抛物型偏微分方程中,初边值问题可以表示为:u(x, 0) = u_0(x)u(x, t) = u_t(x, t) 在边界x=0,x=L上三、求解方法1.紧差分格式紧差分格式是一种求解偏微分方程的方法,其精度为O(h^(1/2) * Δt),无条件稳定。
在这种方法中,我们首先需要建立离散的网格系统,然后通过数值积分求解离散化的偏微分方程。
2.追赶法追赶法是一种求解线性方程组的方法,也可以用于求解初边值问题。
在这种方法中,我们首先需要将偏微分方程转化为线性方程组,然后使用追赶法求解线性方程组。
3.有限元算法有限元算法是一种基于变分原理的求解方法,可以将偏微分方程问题转化为求解有限元空间的线性方程组。
这种方法在求解一维抛物型偏微分方程时具有较高的精度和可靠性。
(完整版)偏微分方程的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程序)
(偏)微分方程一类边值问题的数值求解本文介绍了椭圆型(偏)微分方程一类边值问题的数值求解程序(笔者自编)及其使用方法。
程序基于差分原理,将连续的(偏)微分方程在节点上离散化,最终化为线性代数方程组。
对于原理方面的问题很多微分方程数值解的参考书中都有详尽的描述,但一般缺少具体的实现程序。
笔者认为,对这类数值方法是否理解并能运用的关键之一还是在于是否能写出用于解决问题的程序。
理论基础固然必不可少,程序确往往是我们解决问题的敲门砖。
一维椭圆型方程可表示如下:],[,0)(,0)()(b a x x q x p f qu dxdur dx du p dx d Lu ∈>>=++-=其中L 表示微分算子,很明显这是一个线性算子。
这里要求p ,q 大于零是为了保证最终得到的线性方程组有唯一的非零解,但事实上不满足这个条件可能也是有解的,这涉及到微分方程解的存在性也确定性问题,读者若有兴趣可参考相关书籍。
更具体的,我们可以举出一个一维椭圆型方程的例子来: 例(1):)cos(10)exp(222x u x dx dux dxu d =++- 即:)cos(10),ex p()(,)(,1)(2x f x x q x x r x p ====⎩⎨⎧==3)3(1)1(u u 利用文中提及的程序,可以将以上问题表述为:syms x %定义符号变量xp=@(x) 1; r=@(x) x.^2; q=@(x) exp(x); f=@(x) 10*cos(x);Odbvp(p,r,q,f,1,3,1,3);11.2 1.4 1.6 1.82 2.2 2.4 2.6 2.83-0.50.511.522.532~Odbvp u=f(x)xu哈哈,一条非常漂亮的曲线。
若不满足p,q>0,我们也举一例: 例(2)syms xp=@(x) 10*cos(x); r=@(x) x.^2;q=@(x) 10*sin(x); f=@(x) 10*cos(x);Odbvp(p,r,q,f,1,3,1,3)1 1.2 1.4 1.6 1.82 2.2 2.4 2.6 2.83-500501001502002503003502~Odbvp u=f(x)xu可以发现,程序仍能求解,但结果的光滑性不好。
matlab解偏微分方程
ui,j +1 − ui,j = Hui,j +1 ∆t Hui,j = a2 ui+1,j − 2ui,j + ui−1,j (∆x)2
ui,j +1 − ui,j = Hui,j ∆t 将显式与隐式相加,得平均公式 ui,j +1 − ui,j 1 1 = Hui,j + Hui,j +1 ∆t 2 2
得ui,0 = ui,2 − 2ψi
1 ui,2 = [c(ui+1,1 + ui−1,1) + 2(1 − c)ui,1 + 2ψi t] 2
3.3
例题 两端固定的弦振动
两端固定的弦, 初速为零,初位移是 h x, (0 ≤ x ≤ 2/3) 2 / 3 u(x, 0) = 1−x , (2/3 < x ≤ 1) h 1 − 2/3
作图所用程序如下,其中取c = 0.05, l = 1, h = 0.05.这里使用的方程 与初始条件表示方法与上一节相同. N=4000; c=0.05; x=linspace(0,1,420)’; u1(1:420)=0; u2(1:420)=0; u3(1:420)=0; u1(2:280)=0.05/279*(1:279)’; u1(281:419)=0.05/(419-281)*(419-(281:419)’); u2(2:419)=u1(2:419)+c/2*(u1(3:420)-2*u1(2:419)+u1(1:418)); h=plot(x,u1,’linewidth’,3); axis([0,1,-0.05,0.05]); set(h,’EraseMode’,’xor’,’MarkerSize’,18) for k=2:N set(h,’XData’,x,’YData’,u2) ; drawnow; u3(2:419)=2*u2(2:419)-u1(2:419)+c*(u2(3:420)... -2*u2(2:419)+u2(1:418)); u1=u2; u2=u3; end
【精品】偏微分的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
matlab 偏微分方程
MATLAB是一个强大的数值计算环境,可以用来解决各种各样的数学问题,包括偏微分方程。
下面是一个简单的例子,展示如何在MATLAB中解决一维的偏微分方程。
假设我们要解决以下一维的热传导方程:
∂u∂t=∂2u∂x2
在给定的初始条件和边界条件下:
u(x,0)=sin(πx)u(0,t)=0, u(1,t)=0
我们可以使用MATLAB中的pdepe函数来求解这个问题。
以下是一个简单的MATLAB代码示例:
```matlab
定义参数
T = 1; 最终时间
h = 0.01; 空间步长
t = 0:T/h:T; 时间向量
x = 0:h:1; 空间向量
n = length(x); 空间点的数量
m = length(t); 时间点的数量
初始化矩阵存储解
U = zeros(m, n);
U(:,1) = sin(pi*x); 初始条件
定义偏微分方程
pdepe('u_tt', U, t, x, 'heat', 'periodic');
使用pdepe求解偏微分方程
[U, ~] = pdepe(U, t, x);
绘制结果
surf(x, t, U);
```
这个代码示例使用了MATLAB的pdepe函数,这是一个用于求解偏微分方程的函数。
在上面的代码中,我们首先定义了参数,然后初始化了存储解的矩阵。
然后,我们定义了偏微分方程,并使用pdepe 函数求解它。
最后,我们使用surf函数绘制了结果。
Matlab偏微分方程求解方法
Matlab 偏微分方程求解方法目录:§1 Function Summary on page 10-87§2 Initial Value Problems on page 10-88§3 PDE Solver on page 10-89§4 Integrator Options on page 10-92§5 Examples” on page 10-93§1 Function Summary1.1 PDE Solver” on page 10-871,2 PDE Helper Functi on” on page 10-871.3 PDE SolverThis is the MATLAB PDE solver.PDE Helper Function§2 Initial Value Problemspdepe solves systems of parabolic and elliptic PDEs in one spatial variable x and time t, of the form)xu ,u ,t ,x (s ))x u ,u ,t ,x (f x (x x t u )x u ,u ,t ,x (c m m ∂∂+∂∂∂∂=∂∂∂∂- (10-2) The PDEs hold for b x a ,t t t f 0≤≤≤≤.The interval [a, b] must be finite. mcan be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry,respectively. If m > 0, thena ≥0 must also hold.In Equation 10-2,)x /u ,u ,t ,x (f ∂∂ is a flux term and )x /u ,u ,t ,x (s ∂∂ is a source term. The flux term must depend on x /u ∂∂. The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix )x /u ,u ,t ,x (c ∂∂. The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic equation and otherwise to a parabolic equation. There must be at least one parabolic equation. An element of c that corresponds to a parabolic equation can vanish at isolated values of x if they are mesh points.Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface.At the initial time t = t0, for all x the solution components satisfy initial conditions of the form)x (u )t ,x (u 00= (10-3)At the boundary x = a or x = b, for all t the solution components satisfy a boundary condition of the form0)xu ,u ,t ,x (f )t ,x (q )u ,t ,x (p =∂∂+ (10-4) q(x, t) is a diagonal matrix with elements that are either identically zero or never zero. Note that the boundary conditions are expressed in terms of the f rather than partial derivative of u with respect to x-x /u ∂∂. Also, ofthe two coefficients, only p can depend on u.§3 PDE Solver3.1 The PDE SolverThe MATLAB PDE solver, pdepe, solves initial-boundary value problems for systems of parabolic and elliptic PDEs in the one space variable x and time t.There must be at least one parabolic equation in the system.The pdepe solver converts the PDEs to ODEs using a second-order accurate spatial discretization based on a fixed set of user-specified nodes. The discretization method is described in [9]. The time integration is done with ode15s. The pdepe solver exploits the capabilities of ode15s for solving the differential-algebraic equations that arise when Equation 10-2 contains elliptic equations, and for handling Jacobians with a specified sparsity pattern. ode15s changes both the time step and the formula dynamically.After discretization, elliptic equations give rise to algebraic equations. If the elements of the initial conditions vector that correspond to elliptic equations are not “consistent” with the discretization, pdepe tries to adjust them before eginning the time integration. For this reason, the solution returned for the initial time may have a discretization error comparable to that at any other time. If the mesh is sufficiently fine, pdepe can find consistent initial conditions close to the given ones. If pdepe displays amessage that it has difficulty finding consistent initial conditions, try refining the mesh. No adjustment is necessary for elements of the initial conditions vector that correspond to parabolic equations.PDE Solver SyntaxThe basic syntax of the solver is:sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)Note Correspondences given are to terms used in “Initial Value Problems” on page 10-88.The input arguments arem: Specifies the symmetry of the problem. m can be 0 =slab, 1 = cylindrical, or 2 = spherical. It corresponds to m in Equation 10-2. pdefun: Function that defines the components of the PDE. Itcomputes the terms f,c and s in Equation 10-2, and has the form[c,f,s] = pdefun(x,t,u,dudx)where x and t are scalars, and u and dudx are vectors that approximate the solution and its partial derivative with respect to . c, f, and s are column vectors. c stores the diagonal elements of the matrix .icfun: Function that evaluates the initial conditions. It has the formu = icfun(x)When called with an argument x, icfun evaluates and returns the initial values of the solution components at x in the column vector u.bcfun:Function that evaluates the terms and of the boundary conditions. Ithas the form[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)where ul is the approximate solution at the left boundary xl = a and ur is the approximate solution at the right boundary xr = b. pl and ql are column vectors corresponding to p and the diagonal of q evaluated at xl. Similarly, pr and qr correspond to xr. When m>0 and a = 0, boundedness of the solution near x = 0 requires that the f vanish at a = 0. pdepe imposes this boundary condition automatically and it ignores values returned in pl and ql.xmesh:Vector [x0, x1, ..., xn] specifying the points at which a numerical solution is requested for every value in tspan. x0 and xn correspond to a and b , respectively. Second-order approximation to the solution is made on the mesh specified in xmesh. Generally, it is best to use closely spaced mesh points where the solution changes rapidly. pdepe does not select the mesh in automatically. You must provide an appropriate fixed mesh in xmesh. The cost depends strongly on the length of xmesh. When , it is not necessary to use a fine mesh near to x=0 account for the coordinate singularity.The elements of xmesh must satisfy x0 < x1 < ... < xn.The length of xmesh must be ≥3.tspan:Vector [t0, t1, ..., tf] specifying the points at which a solution is requested for every value in xmesh. t0 and tf correspond tot and f t,respectively.pdepe performs the time integration with an ODE solver that selects both the time step and formula dynamically. The solutions at the points specified in tspan are obtained using the natural continuous extension of the integration formulas. The elements of tspan merely specify where you want answers and the cost depends weakly on the length of tspan.The elements of tspan must satisfy t0 < t1 < ... < tf.The length of tspan must be ≥3.The output argument sol is a three-dimensional array, such that•sol(:,:,k) approximates component k of the solution .•sol(i,:,k) approximates component k of the solution at time tspan(i) and mesh points xmesh(:).•sol(i,j,k) approximates component k of the solution at time tspan(i) and the mesh point xmesh(j).4.2 PDE Solver OptionsFor more advanced applications, you can also specify as input arguments solver options and additional parameters that are passed to the PDE functions.options:Structure of optional parameters that change the default integration properties. This is the seventh input argument.sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)See “Integrator Options” on page 10-92 for more information.Integrator OptionsThe default integration properties in the MATLAB PDE solver are selected to handle common problems. In some cases, you can improve solver performance by overriding these defaults. You do this by supplying pdepe with one or more property values in an options structure.sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)Use odeset to create the options structure. Only those options of the underlying ODE solver shown in the following table are available for pdepe.The defaults obtained by leaving off the input argument options are generally satisfactory. “Integrator Options” on page 10-9 tells you how to create the structure and describes the properties.PDE Properties§4 Examples•“Single PDE” on page 10-93•“System of PDEs” on page 10-98•“Additional Examples” on page 10-1031.Single PDE• “Solving the Equation” on page 10-93• “Evaluating the Solution” on page 10-98Solving the Equation. This example illustrates the straightforward formulation, solution, and plotting of the solution of a single PDE222x u t u ∂∂=∂∂π This equation holds on an interval 1x 0≤≤ for times t ≥ 0. At 0t = the solution satisfies the initial condition x sin )0,x (u π=.At 0x =and 1x = , the solution satisfies the boundary conditions0)t ,1(xu e ,0)t ,0(u t =∂∂+π=- Note The demo pdex1 contains the complete code for this example. The demo uses subfunctions to place all functions it requires in a single MATLAB file.To run the demo type pdex1 at the command line. See “PDE Solver Syntax” on page 10-89 for more information. 1 Rewrite the PDE. Write the PDE in the form)xu ,u ,t ,x (s ))x u ,u ,t ,x (f x (x x t u )x u ,u ,t ,x (c m m ∂∂+∂∂∂∂=∂∂∂∂- This is the form shown in Equation 10-2 and expected by pdepe. For this example, the resulting equation is0x u x x x t u 002+⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂π with parameter and the terms 0m = and the term0s ,xu f ,c 2=∂∂=π=2 Code the PDE. Once you rewrite the PDE in the form shown above (Equation 10-2) and identify the terms, you can code the PDE in a function that pdepe can use. The function must be of the form[c,f,s] = pdefun(x,t,u,dudx)where c, f, and s correspond to the f ,c and s terms. The code below computes c, f, and s for the example problem.function [c,f,s] = pdex1pde(x,t,u,DuDx)c = pi^2;f = DuDx;s = 0;3 Code the initial conditions function. You must code the initial conditions in a function of the formu = icfun(x)The code below represents the initial conditions in the function pdex1ic. Partial Differential Equationsfunction u0 = pdex1ic(x)u0 = sin(pi*x);4 Code the boundary conditions function. You must also code the boundary conditions in a function of the form[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t) The boundary conditions, written in the same form as Equation 10-4, are0x ,0)t ,0(x u .0)t ,0(u ==∂∂+and1x ,0)t ,1(xu .1e t ==∂∂+π- The code below evaluates the components )u ,t ,x (p and )u ,t ,x (q of the boundary conditions in the function pdex1bc.function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)pl = ul;ql = 0;pr = pi * exp(-t);qr = 1;In the function pdex1bc, pl and ql correspond to the left boundary conditions (x=0 ), and pr and qr correspond to the right boundary condition (x=1).5 Select mesh points for the solution. Before you use the MATLAB PDE solver, you need to specify the mesh points at which you want pdepe to evaluate the solution. Specify the points as vectors t and x.The vectors t and x play different roles in the solver (see “PDE Solver” on page 10-89). In particular, the cost and the accuracy of the solution depend strongly on the length of the vector x. However, the computation is much less sensitive to the values in the vector t.10 CalculusThis example requests the solution on the mesh produced by 20 equally spaced points from the spatial interval [0,1] and five values of t from thetime interval [0,2].x = linspace(0,1,20);t = linspace(0,2,5);6 Apply the PDE solver. The example calls pdepe with m = 0, the functions pdex1pde, pdex1ic, and pdex1bc, and the mesh defined by x and t at which pdepe is to evaluate the solution. The pdepe function returns the numerical solution in a three-dimensional array sol, wheresol(i,j,k) approximates the kth component of the solution,u, evaluated atkt(i) and x(j).m = 0;sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);This example uses @ to pass pdex1pde, pdex1ic, and pdex1bc as function handles to pdepe.Note See the function_handle (@), func2str, and str2func reference pages, and the @ section of MATLAB Programming Fundamentals for information about function handles.7 View the results. Complete the example by displaying the results:a Extract and display the first solution component. In this example, the solution has only one component, but for illustrative purposes, the example “extracts” it from the three-dimensional array. The surface plot shows the behavior of the solution.u = sol(:,:,1);surf(x,t,u)title('Numerical solution computed with 20 mesh points') xlabel('Distance x') ylabel('Time t')Distance xNumerical solution computed with 20 mesh points.Time tb Display a solution profile at f t , the final value of . In this example,2t t f ==.figure plot(x,u(end,:)) title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)')Solutions at t = 2.Distance xu (x ,2)Evaluating the Solution. After obtaining and plotting the solution above, you might be interested in a solution profile for a particular value of t, or the time changes of the solution at a particular point x. The kth column u(:,k)(of the solution extracted in step 7) contains the time history of the solution at x(k). The jth row u(j,:) contains the solution profile at t(j). Using the vectors x and u(j,:), and the helper function pdeval, you can evaluate the solution u and its derivative at any set of points xout [uout,DuoutDx] = pdeval(m,x,u(j,:),xout)The example pdex3 uses pdeval to evaluate the derivative of the solution at xout = 0. See pdeval for details.2. System of PDEsThis example illustrates the solution of a system of partial differential equations. The problem is taken from electrodynamics. It has boundary layers at both ends of the interval, and the solution changes rapidly for small . The PDEs are)u u (F xu017.0t u )u u (F xu 024.0t u 212222212121-+∂∂=∂∂--∂∂=∂∂ where )y 46.11exp()y 73.5exp()y (F --=. The equations hold on an interval1x 0≤≤ for times 0t ≥.The solution satisfies the initial conditions0)0,x (u ,1)0,x (u 21≡≡and boundary conditions0)t ,1(xu,0)t ,1(u ,0)t ,0(u ,0)t ,0(x u 2121=∂∂===∂∂ Note The demo pdex4 contains the complete code for this example. The demo uses subfunctions to place all required functions in a single MATLAB file. To run this example type pdex4 at the command line.1 Rewrite the PDE. In the form expected by pdepe, the equations are⎥⎦⎤⎢⎣⎡---+⎥⎦⎤⎢⎣⎡∂∂∂∂∂∂=⎥⎦⎤⎢⎣⎡∂∂⎥⎦⎤⎢⎣⎡)u u (F )u u (F )x /u 170.0)x /u (024.0x u u t *.1121212121 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⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡∂∂∂∂⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡00)x /u (170.0)x /u (024.0*.01u 0212and the right boundary condition is⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡∂∂∂∂⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡-00)x /u (170.0)x /u (024.0*.1001u 2112 Code the PDE. After you rewrite the PDE in the form shown above, you can code it as a function that pdepe can use. The function must be of the form[c,f,s] = pdefun(x,t,u,dudx)where c, f, and s correspond to the , , and terms in Equation 10-2. 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];3 Code the initial conditions function. The initial conditions function must be of the formu = icfun(x)The code below represents the initial conditions in the function pdex4ic. function u0 = pdex4ic(x);u0 = [1; 0];4 Code the boundary conditions function. The boundary conditions functions must be of the form[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)The code below evaluates the components p(x,t,u) and q(x,t) (Equation 10-4) of the boundary conditions in the function pdex4bc.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];5 Select mesh points for the solution. The solution changes rapidly for small t . The program selects the step size in time to resolve this sharp change, but to see this behavior in the plots, output times must be selected accordingly. There are boundary layers in the solution at both ends of [0,1], so mesh points must be placed there to resolve these sharp changes. Often some experimentation is needed to select the mesh that reveals the behavior of the solution.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];6 Apply the PDE solver. The example calls pdepe with m = 0, the functions pdex4pde, pdex4ic, and pdex4bc, and the mesh defined by x and t at which pdepe is to evaluate the solution. The pdepe function returns the numerical solution in a three-dimensional array sol, wheresol(i,j,k) approximates the kth component of the solution, μk, evaluated at t(i) and x(j).m = 0;sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);7 View the results. The surface plots show the behavior of the solution components. u1 = sol(:,:,1); u2 = sol(:,:,2); figure surf(x,t,u1) title('u1(x,t)') xlabel('Distance x') 其输出图形为Distance xu1(x,t)Time tfigure surf(x,t,u2) title('u2(x,t)') xlabel('Distance x') ylabel('Time t')Distance xu2(x,t)Time tAdditional ExamplesThe following additional examples are available. Type edit examplename to view the code and examplename to run the example.。
【VIP专享】一维抛物线偏微分方程数值解法(4)(附图及matlab程序)
一维抛物线偏微分方程数值解法(4)上一篇参看一维抛物线偏微分方程数值解法(3)(附图及matlab程序)解一维抛物线型方程(理论书籍可以参看孙志忠:偏微分方程数值解法)Ut-Uxx=0, 0<x<1,0<t<=1(Ut-aUxx=f(x,t),a>0)U(x,0)=e^x, 0<=x<=1,U(0,t)=e^t,U(1,t)=e^(1+t), 0<t<=1精确解为:U(x,t)=e^(x+t);用紧差分格式:此种方法精度为o(h1^2+h2^4),无条件差分稳定;一:用追赶法解线性方程组(还可以用迭代法解)Matlab程序为:function [u p e x t]=JCHGS(h1,h2,m,n)%紧差分格式解一维抛物线型偏微分方程%此程序用的是追赶法解线性方程组%h1为空间步长,h2为时间步长%m,n分别为空间,时间网格数%p为精确解,u为数值解,e为误差x=(0:m)*h1+0; x0=(0:m)*h1;%定义x0,t0是为了f(x,t)~=0的情况%t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;syms f;for(i=1:n+1)for(j=1:m+1)f(i,j)=0; %f(i,j)=f(x0(j),t0(i))==0%endendfor(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endr=h2/(h1*h1);for(i=1:n) %外循环,先固定每一时间层,每一时间层上解一线性方程组% a(1)=0;b(1)=5/6+r;c(1)=1/12-r/2;d(1)=(r/2-1/12)*u(i+1,1)+...(1/12+r/2)*u(i,1)+(5/6-r)*u(i,2)+(1/12+r/2)*u(i,3)+...h2/12*(f(i,1)+10*f(i,2)+f(i,3));for(k=2:m-2)a(k)=1/12-r/2;b(k)=5/6+r;c(k)=1/12-r/2;d(k)=h2/12*(f(i,k)+...10*f(i,k+1)+f(i,k+2))+(1/12+r/2)*(u(i,k)+u(i,k+2))+(5/6-r)...*u(i,k+1);%输入部分系数矩阵,为0的矩阵元素不输入%一定要注意输入元素的正确性enda(m-1)=1/12-r/2;b(m-1)=5/6+r;d(m-1)=(1/12+r/2)*(u(i,m-1)+u(i,m+1))+...(5/6-r)*u(i,m)+(r/2-1/12)*u(i+1,m+1)+ ...h2/12*(f(i,m-1)+10*f(i,m)+f(i,m+1));for(k=1:m-2) %开始解线性方程组消元过程a(k+1)=-a(k+1)/b(k);b(k+1)=b(k+1)+a(k+1)*c(k);d(k+1)=d(k+1)+a(k+1)*d(k);endu(i+1,m)=d(m-1)/b(m-1); %回代过程%for(k=m-2:-1:1)u(i+1,k+1)=(d(k)-c(k)*u(i+1,k+2))/b(k);endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j)+t(i)); %p为精确解e(i,j)=abs(u(i,j)-p(i,j));%e为误差endend[u p e x t]=JCHGS(0.1,0.005,10,200); surf(x,t,e)>> title('误差');运行约43秒;[u p e x t]=JCHGS(0.1,0.01,10,100);surf(x,t,e) 20多秒;[u p e x t]=JCHGS(0.2,0.04,5,25);surf(x,t,e) 3秒;此方法精度很高;二:g-s迭代法求解线性方程组Matlab程序function [u e p x t k]=JCFGS1(h1,h2,m,n,kmax,ep) % 解抛物线型一维方程格式(Ut-aUxx=f(x,t),a>0)%用g-s(高斯-赛德尔)迭代法解%kmax为最大迭代次数%m,n为x,t方向的网格数,例如(2-0)/0.01=200;%e为误差,p为精确解syms temp;u=zeros(n+1,m+1);x=0+(0:m)*h1;t=0+(0:n)*h2;for(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endfor(i=1:n+1)for(j=1:m+1)f(i,j)=0;endenda=zeros(n,m-1);r=h2/(h1*h1); %此处r=a*h2/(h1*h1);a=1for(k=1:kmax)for(i=1:n)for(j=2:m)temp=((1/12+r/2)*(u(i,j-1)+u(i,j+1))+(5/6-r)*u(i,j)+...h2/12*(f(i,j-1)+10*f(i,j)+f(i,j+1))+(r/2-1/12)*(u(i+1,...j-1)+u(i+1,j+1)))/(5/6+r);a(i+1,j)=(temp-u(i+1,j))*(temp-u(i+1,j));u(i+1,j)=temp;%此处注意是u(i+1,j),,而不是u(i+1,j+1)% endenda(i+1,j)=sqrt(a(i+1,j));if(k>kmax)break;endif(max(max(a))<ep)break;endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j)+t(i));e(i,j)=abs(u(i,j)-p(i,j));endend[u e p x t k]=JCFGS1(0.1,0.005,10,200,100000,1e-12);k=67;运行速度1秒左右;surf(x,t,e)[u e p x t k]=JCFGS1(0.01,0.001,100,1000,1000000,1e-12);k=5780;surf(x,t,e)。
【VIP专享】一维抛物线偏微分方程数值解法(3)(附图及matlab程序)
一维抛物线偏微分方程数值解法(3)上一篇参看一维抛物线偏微分方程数值解法(2)(附图及matlab程序)解一维抛物线型方程(理论书籍可以参看孙志忠:偏微分方程数值解法)Ut-Uxx=0, 0<x<1,0<t<=1(Ut-aUxx=f(x,t),a>0)U(x,0)=e^x, 0<=x<=1,U(0,t)=e^t,U(1,t)=e^(1+t), 0<t<=1精确解为:U(x,t)=e^(x+t);此种方法精度为o(h1^2+h2^2)一:用追赶法解线性方程组(还可以用迭代法解)Matlab程序function [u p e x t]=CN(h1,h2,m,n)%Crank-Nicolson格式差分法解一维抛物线型偏微分方程%此程序用的是追赶法解线性方程组%h1为空间步长,h2为时间步长%m,n分别为空间,时间网格数%p为精确解,u为数值解,e为误差x=(0:m)*h1+0; x0=(0:m)*h1; %定义x0,t0是为了f(x,t)~=0的情况%t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;syms f;for(i=1:n+1)for(j=1:m+1)f(i,j)=0; %f(i,j)=f(x0(j),t0(i))==0%endendfor(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endr=h2/(h1*h1);for(i=1:n) %外循环,先固定每一时间层,每一时间层上解一线性方程组% a(1)=0;b(1)=1+r;c(1)=-r/2;d(1)=r/2*(u(i+1,1)+u(i,1))+h2*f(i,j)...+(1-r)*u(i,2)+r/2*u(i,3);for(k=2:m-2)a(k)=-r/2;b(k)=1+r;c(k)=-r/2;d(k)=h2*f(i,j)+r/2*u(i,k)+(1-r)...*u(i,k+1)+r/2*u(i,k+2);%输入部分系数矩阵,为0的矩阵元素不输入%enda(m-1)=-r/2;b(m-1)=1+r;d(m-1)=h2*f(i,j)+r/2*(u(i,m+1)+u(i+1,m+1)...)+r/2*u(i,m-1)+(1-r)*u(i,m);for(k=1:m-2) %开始解线性方程组消元过程a(k+1)=-a(k+1)/b(k);b(k+1)=b(k+1)+a(k+1)*c(k);d(k+1)=d(k+1)+a(k+1)*d(k);endu(i+1,m)=d(m-1)/b(m-1); %回代过程%for(k=m-2:-1:1)u(i+1,k+1)=(d(k)-c(k)*u(i+1,k+2))/b(k);endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j)+t(i)); %p为精确解e(i,j)=abs(u(i,j)-p(i,j));%e为误差endend[u p e x t]=CN(0.1,0.005,10,200);surf(x,t,e); shading interp;>> xlabel('x');ylabel('t');zlabel('e');>> title('误差曲面')plot(x,e)plot(t,e)误差较向前欧拉法减小一半但是运行时间较长,约39秒,而前两次运行只需l秒左右;[u p e x t]=CN(0.01,0.01,100,100);运行需三分钟左右,误差比前次提高五倍,运算量也提高五倍[u p e x t]=CN(0.1,0.1,10,10);surf(x,t,e) 运行需要2秒;精度还是挺高的;[u p e x t]=CN(0.1,0.2,10,5);surf(x,t,e)误差还可以接受此种方法精度高,计算量较大二:用迭代法解线性方程组:Matlab程序如下:function [u e p x t k]=CN1(h1,h2,m,n,kmax,ep) % 解抛物线型一维方程 C-N格式(Ut-aUxx=f(x,t),a>0) %用g-s(高斯-赛德尔)迭代法解%kmax为最大迭代次数%m,n为x,t方向的网格数,例如(2-0)/0.01=200;%e为误差,p为精确解syms temp;u=zeros(n+1,m+1);x=0+(0:m)*h1;t=0+(0:n)*h2;for(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endfor(i=1:n+1)for(j=1:m+1)f(i,j)=0;endenda=zeros(n,m-1);r=h2/(h1*h1); %此处r=a*h2/(h1*h1);a=1for(k=1:kmax)for(i=1:n)for(j=2:m)temp=((r/2*u(i,j-1)+(1-r)*u(i,j)+r/2*u(i,...j+1)+h2*f(i,j)+r/2*u(i+1,j-1)+r/2*u(i+1,j+1))/(1+r));a(i+1,j)=(temp-u(i+1,j))*(temp-u(i+1,j));u(i+1,j)=temp;%此处注意是u(i+1,j),,而不是u(i+1,j+1)% endenda(i+1,j)=sqrt(a(i+1,j));if(k>kmax)break;endif(max(max(a))<ep)break;endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j)+t(i));e(i,j)=abs(u(i,j)-p(i,j));endend[u e p x t k]=CN1(0.1,0.005,10,200,10000,1e-10);运行速度:1秒迭代次数k =81surf(x,t,e)第二幅图为三角追赶法解方程作出的图,两者几乎一样;由于迭代法速度很快,所以可以将区间分得更小[u e p x t k]=CN1(0.01,0.01,100,100,10000,1e-12);surf(x,t,e);shading interp; k=6903。
偏微分方程数值解法的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中的偏微分方程数值解法
MATLAB中的偏微分方程数值解法偏微分方程(Partial Differential Equations,PDEs)是数学中的重要概念,广泛应用于物理学、工程学、经济学等领域。
解决偏微分方程的精确解往往非常困难,因此数值方法成为求解这类问题的有效途径。
而在MATLAB中,有丰富的数值解法可供选择。
本文将介绍MATLAB中几种常见的偏微分方程数值解法,并通过具体案例加深对其应用的理解。
一、有限差分法(Finite Difference Method)有限差分法是最为经典和常用的偏微分方程数值解法之一。
它将偏微分方程的导数转化为差分方程,通过离散化空间和时间上的变量,将连续问题转化为离散问题。
在MATLAB中,使用有限差分法可以比较容易地实现对偏微分方程的数值求解。
例如,考虑一维热传导方程(Heat Equation):∂u/∂t = k * ∂²u/∂x²其中,u为温度分布随时间和空间的变化,k为热传导系数。
假设初始条件为一段长度为L的棒子上的温度分布,边界条件可以是固定温度、热交换等。
有限差分法可以将空间离散化为N个节点,时间离散化为M个时刻。
我们可以使用中心差分近似来计算二阶空间导数,从而得到以下差分方程:u(i,j+1) = u(i,j) + Δt * (k * (u(i+1,j) - 2 * u(i,j) + u(i-1,j))/Δx²)其中,i表示空间节点,j表示时间步。
Δt和Δx分别为时间和空间步长。
通过逐步迭代更新节点的温度值,我们可以得到整个时间范围内的温度分布。
而MATLAB提供的矩阵计算功能,可以大大简化有限差分法的实现过程。
二、有限元法(Finite Element Method)有限元法是另一种常用的偏微分方程数值解法,特点是适用于复杂的几何形状和边界条件。
它将求解区域离散化为多个小单元,通过构建并求解代数方程组来逼近连续问题。
在MATLAB中,我们可以使用Partial Differential Equation Toolbox提供的函数进行有限元法求解。
偏微分方程数值解法(3)
§4 双曲型方程的差分解法 一、一阶双曲型方程的差分格式 一阶双曲型方程的初值问题为⎩⎨⎧=>+∞<<-∞=+)2()()0,()1()0,(0x x u t x au u x t ϕa 为常数,亦称(1)为对流方程。
称ξ=-at x 为(1)的特征线,ξ 为常数,沿特征线 u (x , t )的方向导数0d ),(d d =+=+=t x u au tt at du t u ξ 即u (x , t )沿特征线为常数,再由 (2.),得初值问题(1),(2)的解)(),(at x t x u -=ϕ这是个单向的传播波,a >0时,波形ϕ(x )沿x 轴方向传播,为右传播波,a < 0时,为左传播波,在传播过程中,波形均不发生变化。
二阶波动方程02=-xx tt u a u若令v = u ,w = au x 则得一阶双曲型方程组⎩⎨⎧=-=-0x t x t av w aw v 再令v w vv w u-=+=~,~,则得⎩⎨⎧=+=-0~~0~~x tx t v a v u a u 可见,二阶双曲型方程可化为一阶双曲型方程组。
下面建立(1)的差分格式,作网格线 ,2,1,0,±±===j jh x x j,2,1,0,===n n t t n τ对区域G :}0,),({>+∞<<∞-t x t x 进行剖分,其中h = ∆x 为空间步长,τ = ∆t 为时间步长。
a )逆风格式 u t (x j , t n )用向前差商代替,u x (x j , t n )用向前或向后差商代替,nj u 表示u (x j , t n )近似值,得011=-+-++hu u au u njn j n jn j τ或011=-+--+hu u au u nj n j n jn j τ令λ = τ / h ,得 )(11nj n j n j n j u u a u u --=++λ (3))(11n j n j n j n j u u a u u -+--=λ(4)截断误差均为)(h o +τ,其节点分布见图1。
偏微分方程数值解法的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 ~= mdisp('输入参数有误!')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迭代格式。
一维抛物线偏微分方程数值解法
一维抛物线偏微分方程数值解法(2)上一篇文章请参看一维抛物线偏微分方程数值解法(1)解一维抛物线型方程(理论书籍可以参看孙志忠:偏微分方程数值解法)Ut-Uxx=0, 0<x<1,0<t<=1(Ut-aUxx=f(x,t),a>0)U(x,0)=e^x, 0<=x<=1,U(0,t)=e^t,U(1,t)=e^(1+t), 0<t<=1精确解为:U(x,t)=e^(x+t);Matlab程序:(此为向后差分法)function [u p e x t]=pwxywxh(h1,h2,m,n)%欧拉向后差分法解一维抛物线型偏微分方程%此程序用的是追赶法解线性方程组%h1为空间步长,h2为时间步长%m,n分别为空间,时间网格数%p为精确解,u为数值解,e为误差x=(0:m)*h1+0;t=(0:n)*h2+0;for(i=1:n+1)for(j=1:m+1)f(i,j)=0;endendfor(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endr=h2/(h1*h1);for(i=2:n+1) %外循环,先固定每一时间层,每一时间层上解一线性方程组%a(1)=0;b(1)=1+2*r;c(1)=-r;d(1)=u(i-1,2)+h2*f(i,2)+r*u(i,1);for(k=2:m-2)a(k)=-r;b(k)=1+2*r;c(k)=-r;d(k)=u(i-1,k+1)+h2*f(i,k+1);%输入部分系数矩阵,为0的矩阵元素不输入%enda(m-1)=-r;b(m-1)=1+2*r;d(m-1)=u(i-1,m)+h2*f(i,m)+r*u(i,m+1);for(k=1:m-2) %开始解线性方程组消元过程a(k+1)=-a(k+1)/b(k);b(k+1)=b(k+1)+a(k+1)*c(k);d(k+1)=d(k+1)+a(k+1)*d(k);endu(i,m)=d(m-1)/b(m-1); %回代过程%for(k=m-2:-1:1)u(i,k+1)=(d(k)-c(k)*u(i,k+2))/b(k);endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j)+t(i)); %p为精确解e(i,j)=abs(u(i,j)-p(i,j));%e为误差endend[u p e x t]=pwxywxh(0.1,0.005,10,200);surf(x,t,e);xlabel('x');ylabel('t');zlabel('e');>> title('误差曲面');plot(t,e)误差较之前的欧拉向前差分格式增长了两倍[u p e x t]=pwxywxh(0.1,0.05,10,20); plot(t,e)[u p e x t]=pwxywxh(0.01,0.05,100,20); plot(t,e)[u p e x t]=pwxywxh(0.01,0.005,100,200);plot(x,e)[u p e x t]=pwxywxh(0.005,0.005,200,200); plot(x,e)X=1时,出现了误差??? 不是边界条件吗?不能理解这方法还是比前一种方法误差大呀不过可以随便改变时间、空间步长。
偏微分方程的数值解方法及源程序
-240-第二十章 偏微分方程的数值解自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。
这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。
我们将只含有未知多元函数及其偏导数的方程,称之为偏微分方程。
方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。
如果方程中对于未知函数和它的所有偏导数都是线性的,这样的方程称为线性偏微分方程,否则称它为非线性偏微分方程。
初始条件和边界条件称为定解条件,未附加定解条件的偏微分方程称为泛定方程。
对于一个具体的问题,定解条件与泛定方程总是同时提出。
定解条件与泛定方程作为一个整体,称为定解问题。
§1 偏微分方程的定解问题各种物理性质的定常(即不随时间变化)过程,都可用椭圆型方程来描述。
其最典型、最简单的形式是泊松(Poisson)方程),(2222y x f y ux u u =∂∂+∂∂=Δ (1)特别地,当0),(≡y x f 时,即为拉普拉斯(Laplace)方程,又称为调和方程02222=∂∂+∂∂=Δyux u u (2)带有稳定热源或内部无热源的稳定温度场的温度分布,不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。
Poisson 方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(|),(),(),(),(2222y x y x u y x y x f y uxu y x ϕ (3)其中Ω为以Γ为边界的有界区域,Γ为分段光滑曲线,ΓΩU 称为定解区域,),(),,(y x y x f ϕ分别为ΓΩ,上的已知连续函数。
第二类和第三类边界条件可统一表示成),(),(y x u n u y x ϕα=⎟⎠⎞⎜⎝⎛+∂∂Γ∈ (4) 其中n 为边界Γ的外法线方向。
当0=α时为第二类边界条件,0≠α时为第三类边界条件。
在研究热传导过程,气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时,常常会遇到抛物型方程。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
一维抛物线偏微分方程数值解法(3)
上一篇参看一维抛物线偏微分方程数值解法(2)(附图及matlab程序)
解一维抛物线型方程(理论书籍可以参看孙志忠:偏微分方程数值解法)
Ut-Uxx=0, 0<x<1,0<t<=1(Ut-aUxx=f(x,t),a>0)
U(x,0)=e^x, 0<=x<=1,
U(0,t)=e^t,U(1,t)=e^(1+t), 0<t<=1
精确解为:U(x,t)=e^(x+t);
此种方法精度为o(h1^2+h2^2)
一:用追赶法解线性方程组(还可以用迭代法解)
Matlab程序
function [u p e x t]=CN(h1,h2,m,n)
%Crank-Nicolson格式差分法解一维抛物线型偏微分方程
%此程序用的是追赶法解线性方程组
%h1为空间步长,h2为时间步长
%m,n分别为空间,时间网格数
%p为精确解,u为数值解,e为误差
x=(0:m)*h1+0; x0=(0:m)*h1; %定义x0,t0是为了f(x,t)~=0的情况%
t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;
syms f;
for(i=1:n+1)
for(j=1:m+1)
f(i,j)=0; %f(i,j)=f(x0(j),t0(i))==0% end
end
for(i=1:n+1)
u(i,1)=exp(t(i));
u(i,m+1)=exp(1+t(i));
end
for(i=1:m+1)
u(1,i)=exp(x(i));
end
r=h2/(h1*h1);
for(i=1:n) %外循环,先固定每一时间层,每一时间层上解一线性方程组%
a(1)=0;b(1)=1+r;c(1)=-r/2;d(1)=r/2*
(u(i+1,1)+u(i,1))+h2*f(i,j)...
+(1-r)*u(i,2)+r/2*u(i,3);
for(k=2:m-2)
a(k)=-r/2;b(k)=1+r;c(k)=-
r/2;d(k)=h2*f(i,j)+r/2*u(i,k)+(1-r)...
*u(i,k+1)+r/2*u(i,k+2);
%输入部分系数矩阵,为0的矩阵元素不输入%
end
a(m-1)=-r/2;b(m-1)=1+r;d(m-1)=h2*f(i,j)+r/2*
(u(i,m+1)+u(i+1,m+1)...
)+r/2*u(i,m-1)+(1-r)*u(i,m);
for(k=1:m-2) %开始解线性方程组消元过程
a(k+1)=-a(k+1)/b(k);
b(k+1)=b(k+1)+a(k+1)*c(k);
d(k+1)=d(k+1)+a(k+1)*d(k);
end
u(i+1,m)=d(m-1)/b(m-1); %回代过程%
for(k=m-2:-1:1)
u(i+1,k+1)=(d(k)-c(k)*u(i+1,k+2))/b(k);
end
end
for(i=1:n+1)
for(j=1:m+1)
p(i,j)=exp(x(j)+t(i)); %p为精确解
e(i,j)=abs(u(i,j)-p(i,j));%e为误差
end
end
[u p e x t]=CN(0.1,0.005,10,200);surf(x,t,e); shading interp;
>> xlabel('x');ylabel('t');zlabel('e');
>> title('误差曲面')
plot(x,e)
plot(t,e)
误差较向前欧拉法减小一半
但是运行时间较长,约39秒,而前两次运行只需l秒左右;
[u p e x t]=CN(0.01,0.01,100,100);运行需三分钟左右,误差比前次提高五倍,运算量也提高五倍
[u p e x t]=CN(0.1,0.1,10,10);surf(x,t,e) 运行需要2秒;精度还是挺高的;
[u p e x t]=CN(0.1,0.2,10,5);surf(x,t,e)
误差还可以接受
此种方法精度高,计算量较大
二:用迭代法解线性方程组:
Matlab程序如下:
function [u e p x t k]=CN1(h1,h2,m,n,kmax,ep) % 解抛物线型一维方程 C-N格式(Ut-aUxx=f(x,t),a>0) %用g-s(高斯-赛德尔)迭代法解
%kmax为最大迭代次数
%m,n为x,t方向的网格数,例如(2-0)/0.01=200;
%e为误差,p为精确解
syms temp;
u=zeros(n+1,m+1);
x=0+(0:m)*h1;
t=0+(0:n)*h2;
for(i=1:n+1)
u(i,1)=exp(t(i));
u(i,m+1)=exp(1+t(i));
end
for(i=1:m+1)
u(1,i)=exp(x(i));
end
for(i=1:n+1)
for(j=1:m+1)
f(i,j)=0;
end
end
a=zeros(n,m-1);
r=h2/(h1*h1); %此处r=a*h2/(h1*h1);a=1
for(k=1:kmax)
for(i=1:n)
for(j=2:m)
temp=((r/2*u(i,j-1)+(1-
r)*u(i,j)+r/2*u(i,...
j+1)+h2*f(i,j)+r/2*u(i+1,j-
1)+r/2*u(i+1,j+1))/(1+r));
a(i+1,j)=(temp-u(i+1,j))*(temp-
u(i+1,j));
u(i+1,j)=temp;%此处注意是u(i+1,j),,而不是u(i+1,j+1)%
end
end
a(i+1,j)=sqrt(a(i+1,j));
if(k>kmax)
break;
end
if(max(max(a))<ep)
break;
end
end
for(i=1:n+1)
for(j=1:m+1)
p(i,j)=exp(x(j)+t(i));
e(i,j)=abs(u(i,j)-p(i,j));
end
end
[u e p x t k]=CN1(0.1,0.005,10,200,10000,1e-10);运行速度:1秒迭代次数k =
81
surf(x,t,e)
第二幅图为三角追赶法解方程作出的图,两者几乎一样;
由于迭代法速度很快,所以可以将区间分得更小
[u e p x t k]=CN1(0.01,0.01,100,100,10000,1e-12);surf(x,t,e);shading interp; k=6903。