有限差分求解偏微分方程matlab

如何使用matlab，用有限差分求裂扮解偏微分方程？

求解思路：把偏微分方程离散化，采用合适的差分方法，将复杂的方程简化成简单的线性方程组，最后求解线性方程组，得到其数值解。

现以一维扩散方程为例，说明其计算过程。

第一步，根据条件，建立边界条件和初始条件，即

g0=@(t)zeros(size(t))

g1=g0%边界条件

eta=@(x)sin(pi*x)%初始条件

第二步，设置网格数，即

n=101%网格数

m=101%网格数

第三步，设置步长，即

h=0.01%步长

k=0.01%步长

第四步，设置t和x的初始值，即

t0=0%t的初始值

x0=0%x的初链启始值

第五步，确定扩散系数，即

K=1/pi^2

第六步，自定义Crank-Nicolson差分格式解函数

[t,x,U]=diffusion_sol1(h,k,t0,x0,n,m,eta,g0,g1,K)

第七步，绘制棚源如偏微分方程解的曲面，即

surf(t,x,U)

最后，运行程序得到一维扩散方程数值解的曲面图

用Matlab求解偏微分方程，可以用pde工具箱来解决。例如，简单的点热源方程

求解步骤及主要函数：

1、问题定义

2、创建具有单个腊竖毕因变量的PDE模型纤液， createpde（）

3、创建几何结构并将其追加到PDE模型中，geometryFromEdges（）

4、使用边界条件，pdegplot（）

5、指定PDE系数

6、指定初始条件，setInitialConditions（）

7、生成网格，pdemesh（）

8、生成时间离散化

9、求解数值解，solvepde（）

10、绘图解的图形，pdeplot（）

执行代码：

clc

R1 = [34-111-1-1-111]

C1 = [1000.4]

C1 = [C1zeros(length(R1) - length(C1),1)]

gd = [R1,C1]

sf = 'R1+C1'

ns = char('R1','C1')'

g = decsg(gd,sf,ns)

numberOfPDE = 1

pdem = createpde(numberOfPDE)

geometryFromEdges(pdem,g)

figure

pdegplot(pdem,'edgeLabels','on','subdomainLabels','on')

axis([-1.1 1.1 -1.1 1.1])

axis equal

title 'Geometry With Edge and Subdomain Labels'

applyBoundaryCondition(pdem,'Edge',(1:4),'u',0)

specifyCoefficients(pdem,'m',0, '轮芹d',1, 'c',1, 'a',0,'f',1)

setInitialConditions(pdem,0)

setInitialConditions(pdem,1,'face',2)

msh = generateMesh(pdem)

figure

pdemesh(pdem)

axis equal

nframes = 20

tlist = linspace(0,0.1,nframes)

pdem.SolverOptions.ReportStatistics ='on'

result = solvepde(pdem,tlist)

u1 = result.NodalSolution

figure

umax = max(max(u1))

umin = min(min(u1))

for j = 1:nframes,

pdeplot(pdem,'xydata',u1(:,j),'zdata',u1(:,j))

caxis([umin umax])

axis([-1 1 -1 1 0 1])

Mv(j) = getframe

end

movie(Mv,1)

运行结果

MATLAB提供两种解决PDE问题:

pdepe()函数求解般PDEs据用较通用性支持命令行形式调用

二PDE工具箱求解特殊PDE问题PDEtool较局限性比能求解二阶PDE问题并且能解决偏微程组提供GUI界面繁杂编程解脱同通File->Save As直接M代码

MATLAB语言提供pdepe()函数直接求解般偏微程(组)调用格式

sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)

【输入并渣绝参数】

@pdefun:PDE问题描述函数必须换面标准形式

PDE编写面入口函数

[c,f,s]=pdefun(x,t,u,du)

m,x,t应于(式1)相关参数duu阶导数由给定输入变量即表示c,f,s三函数

@pdebc:PDE边界条件描述函数必梁渗须先化面形式

于边值条件编写面函数描述

[pa,qa,pb,qb]=pdebc(x,t,u,du)

其a表示边界b表示边界

@pdeic:PDE初值条件必须化面形式

股我使用面简单函数描述

u0=pdeic(x)

m,x,t:应于(式1)相关参数

【输参绝姿数】

sol:三维数组sol(:,:,i)表示ui解换句说uk应x(i)t(j)解sol(i,j,k)

通sol我使用pdeval()直接计算某点函数值

---