一维导热方程有限差分法matlab实现

第五次作业（前三题写在作业纸上）

一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf 文件，热扩散系数α=const ，

22T T t x

α??=?? 1. 用Tylaor 展开法推导出FTCS 格式的差分方程

2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。

3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。

4. 编写M 文件求解上述方程，并用适当的文字对程序做出说明。（部分由网络搜索得到，添加，修改后得到。）

function rechuandaopde

%以下所用数据，除了t 的范围我根据题目要求取到了20000，其余均从pdf 中得来 a=0.00001;%a 的取值

xspan=[0 1];%x 的取值范围

tspan=[0 20000];%t 的取值范围

ngrid=[100 10];%分割的份数，前面的是t 轴的，后面的是x 轴的

f=@(x)0;%初值

g1=@(t)100;%边界条件一

g2=@(t)100;%边界条件二

[T,x,t]=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);%计算所调用的函数

[x,t]=meshgrid(x,t);

mesh(x,t,T);%画图，并且把坐标轴名称改为x ，t ，T

xlabel('x')

ylabel('t')

zlabel('T')

T%输出温度矩阵

dt=tspan(2)/ngrid(1);%t 步长

h3000=3000/dt;

h9000=9000/dt;

h15000=15000/dt;%3000,9000,15000下，温度分别在T矩阵的哪些行T3000=T(h3000,:)

T9000=T(h9000,:)

T15000=T(h15000,:)%输出三个时间下的温度分布

%不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面

%稳定性讨论,傅里叶级数法

dx=xspan(2)/ngrid(2);%x步长

sta=4*a*dt/(dx^2)*(sin(pi/2))^2;

if sta>0,sta<2

fprintf('\n%s\n','有稳定性')

else

fprintf('\n%s\n','没有稳定性')

error

end

%真实值计算

[xe,te,Te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid);

[xe,te]=meshgrid(xe,te);

mesh(xe,te,Te);%画图，并且把坐标轴名称改为xe，te，Te

xlabel('xe')

ylabel('te')

zlabel('Te')

Te%输出温度矩阵

%误差计算

jmax=1/dx+1;%网格点数

[rms]=wuchajisuan(T,Te,jmax)

rms%输出误差

function [rms]=wuchajisuan(T,Te,jmax)

for j=1:jmax

rms=((T(j)-Te(j))^2/jmax)^(1/2)

end

function[Ue,xe,te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid)

n=ngrid(1);%t份数

m=ngrid(2);%x份数

Ue=zeros(ngrid);

xe=linspace(xspan(1),xspan(2),m);%画网格

te=linspace(tspan(1),tspan(2),n);%画网格

for j=2:n

for i=2:m-1

for g=1:m-1

Ue(j,i)=100-(400/(2*g-1)/pi)*sin((2*g-1)*pi*xe(j))*exp(-a*(2*g-1)^2*pi^2*te(i)) end

end

end

function [U,x,t]=pdesolution(a,f,g1,g2,xspan,tspan,ngrid)

n=ngrid(1);%t份数

m=ngrid(2);%x份数

h=range(xspan)/(m-1);%x网格长度

x=linspace(xspan(1),xspan(2),m);%画网格

k=range(tspan)/(n-1); %t网格长度

t=linspace(tspan(1),tspan(2),n);%画网格

U=zeros(ngrid);

U(:,1)=g1(t);%边界条件

U(:,m)=g2(t);

U(1,:)=f(x);%初值条件

%差分计算

for j=2:n

for i=2:m -1

U(j,i)=(1-2*a*k/h^2)*U(j -1,i)+a*k/h^2*U(j -1,i -1)+a*k/h^2*U(j -1,i+1);

end

end

5. 将温度随时间变化情况用曲线表示

x t T

6. 给出3000、9000、15000三个时刻的温度分布情况，对温度随时间变化规律做说明。 T3000=100.0000 63.4362 34.2299 15.8021

7.4641 7.4641 15.8021 34.2299 63.4362 100.0000

T9000=100.0000 81.6930 65.6076 53.6839 47.3466 47.3466 53.6839 65.6076 81.6930 100.0000

T15000=100.0000 89.9415 81.0962 74.5310 71.0378 71.0378 74.5310 81.0962 89.9415 100.0000

根据数据分析，在同一个x点上温度随时间的增加而增加，但增幅变小。x-T图形仍为抛物线，但随着时间的增加，极值变小，图像变得平缓。

7.用计算数据说明，并结合差分方程余项，空间、时间间隔对求解精度影响。

数据量较大，且原理相同，我取一个向量演示一下。

保持空间间隔不变，修改时间间隔，时间间隔加大，得到的误差加大。

保持时间间隔不变，修改空间间隔，空间间隔加大，得到的误差加大。

修改空间间隔的误差在增量比修改时间间隔的大。

从方差余项上来看，（没有公式编辑器。。。。只能从ppt里粘贴了）这个余项里的△t，△x都在分母上，所以与误差成正比，且△x的次数应该是比△t高，故影响较大。

8.用计算数据说明，稳定性要求对求解精度的影响。

修改稳定性，即修改x和t分的份数（ngrid），之后看误差。

稳定性越高，解的精度越高。即在满足稳定性要求（a*△t/(△x^2)<0.5）时，a*△t/(△x^2)越接近0，误差越小。

从概念上理解，稳定性越好，对引入时间层误差的抑制能力越强。所以误差越小。

二、调用MATLAB函数完成上述计算

1.编写M文件求解上述方程，并用适当的文字对程序做出说明。

function pdepediaoyong

m=0;

x=linspace(0,1,11);%x的网格

t=linspace(0,20000,101);%t的网格

sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);%调用函数

T=sol(:,:,1);%解

figure;%画图

surf(x,t,T)

xlabel('x')

ylabel('t')

zlabel('T')

dt=20000/100;%t步长

h3000=3000/dt;

h9000=9000/dt;

h15000=15000/dt;%3000,9000,15000下，温度分别在T矩阵的哪些行T3000=T(h3000,:)

T9000=T(h9000,:)

T15000=T(h15000,:)%输出三个时间下的温度分布

%不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面

function [c,f,s]=pdefun(x,t,T,DuDx)%PDE方程函数

c=100000;

f=DuDx;

s=0;

function u0=icfun(x)%初始条件函数

u0=0;

function [pl,ql,pr,qr]=bcfun(xl,Tl,xr,Tr,t)%边界条件函数

pl = Tl-100;

ql = 0;

pr = Tr-100;

qr = 0;

2.将温度随时间变化情况用曲线表示。

3.给出3000、9000、15000三个时刻的温度分布情况。

T3000 =100.0000 67.1058 39.8611 21.1973 10.9885 7.8279 10.9885 21.1973 39.8611 67.1058 100.0000

T9000=100.0000 83.4839 68.6032 56.8191 49.2705 46.6732 49.2705 56.8191 68.6032 83.4839 100.0000

T15000=100.0000 90.8310 82.5601 75.9972 71.7845 70.3330 71.7845 75.9972 82.5601 90.8310 100.0000

根据数据分析，在同一个x点上温度随时间的增加而增加，但增幅变小。x-T图形仍为抛物线，但随着时间的增加，极值变小，图像变得平缓。

4.用计算数据说明，空间、时间间隔对求解精度影响，并与有限差分法的计算结果做

比较。

调用前面做出来的真实值,跟pdepe做出来的值计算误差，再与有限差分法的误差比较。用pdepe函数求的误差远小于有限差分法，所以pdepe函数法更精确。

5.用计算数据说明，有无稳定性要求，为什么？若有，如何对求解精度的影响。

不知道这个pdepe函数的稳定性要用什么检验。傅里叶级数检验不适用。

有限差分法求解偏微分方程MATLAB教学教材

有限差分法求解偏微分方程M A T L A B

南京理工大学 课程考核论文 课程名称：高等数值分析 论文题目：有限差分法求解偏微分方程姓名：罗晨 学号： 115104000545 成绩： 有限差分法求解偏微分方程

一、主要内容 1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程： 22(,)()u u f x t t x αα??-=??其中为常数 具体求解的偏微分方程如下： 22001 (,0)sin()(0,)(1,)00 u u x t x u x x u t u t t π???-=≤≤?????? =??? ==≥??? 2.推导五种差分格式、截断误差并分析其稳定性； 3.编写MATLAB 程序实现五种差分格式对偏微分方程的求解及误差分析； 4.结论及完成本次实验报告的感想。 二、推导几种差分格式的过程： 有限差分法（finite-difference methods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。 推导差分方程的过程中需要用到的泰勒展开公式如下： ()2 100000000()()()()()()()......()(()) 1!2!! n n n f x f x f x f x f x x x x x x x o x x n +'''=+-+-++-+- (2-1) 求解区域的网格划分步长参数如下：

(完整版)大连理工大学高等数值分析抛物型方程有限差分法

抛物型方程有限差分法 1. 简单差分法 考虑一维模型热传导方程 (1.1) )(22x f x u a t u +??=??，T t ≤<0 其中a 为常数。)(x f 是给定的连续函数。（1.1）的定解问题分两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： （1.2） ()()x x u ?=0,, ∞<<∞-x 第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： ()13.1 ()()x x u ?=0,, l x l <<- 及边值条件 ()23.1 ()()0,,0==t l u t u ， T t ≤≤0 假定()x f 和()x ?在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。

现在考虑边值问题（1.1），（1.3）的差分逼近 取 N l h = 为空间步长，M T = τ为时间步长，其中N ，M 是 自然数， jh x x j ==, ()N j ,,1,0Λ=； τ k y y k ==, ()M k ,,1,0Λ= 将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 ()j i y x ,表 示网格节点； h G 表示网格内点（位于开矩形G 中的网格节点）的集合； h G 表示位于闭矩形G 中的网格节点的集合； h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点()k i t x ,处的待求近似解，N j ≤≤0，M k ≤≤0。 注意到在节点()k i t x ,处的微商和差商之间的下列关系 （(,)k j k j u u x t t t ????≡ ? ????）： ()() ()ττ O t u t x u t x u k j k j k j +??? ????=-+,,1 ()() ()2112,,ττ O t u t x u t x u k j k j k j +??? ????=--+ ()()()h O x u h t x u t x u k j k j k j +??? ????=-+,,1 ()() ()h O x u h t x u t x u k j k j k j +??? ????=--,,1 ()() ()2112,,h O x u h t x u t x u k j k j k j +??? ????=--+ ()()() ()2 222 11,,2,h O x u h t x u t x u t x u k j k j k j k j +???? ????=+--+ 可得到以下几种最简差分格式

【毕业设计(论文)】二维热传导方程有限差分法的MATLAB实现

第1章前言 1.1问题背景 在史策教授的《一维热传导方程有限差分法的MATLAB实现》和曹刚教授的《一维偏微分方程的基本解》中，对偏微分方程的解得MATLAB实现问题进行过研究，但只停留在一维中，而实际中二维和三维的应用更加广泛。诸如粒子扩散或神经细胞的动作电位。也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与Ornstein-uhlenbeck过程。热方程及其非线性的推广形式也被应用与影响分析。 在科学和技术发展过程中，科学的理论和科学的实验一直是两种重要的科学方法和手段。虽然这两种科学方法都有十分重要的作用，但是一些研究对象往往由于他们的特性（例如太大或太小，太快或太慢）不能精确的用理论描述或用实验手段来实现。自从计算机出现和发展以来，模拟那些不容易观察到的现象，得到实际应用所需要的数值结果，解释各种现象的规律和基本性质。 科学计算在各门自然科学和技术科学与工程科学中其越来越大的作用，在很多重要领域中成为不可缺少的重要工具。而科学与工程计算中最重要的内容就是求解科学研究和工程技术中出现的各种各样的偏微分方程或方程组。 解偏微分方程已经成为科学与工程计算的核心内容，包括一些大型的计算和很多已经成为常规的计算。为什么它在当代能发挥这样大的作用呢？第一是计算机本身有了很大的发展；第二是数值求解方程的计算法有了很大的发展，这两者对人们计算能力的发展都是十分重要的。 1.2问题现状 近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但在计算数学专业，而且也在其他理工科专业的研究生的大学生中开设。同时，求解热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。而且精度上更好。 目前，在欧美各国MATLAB的使用十分普及。在大学的数学、工程和科学系科，MATLAB

【文献综述】热传导方程差分格式的收敛性和稳定性

文献综述 信息与计算科学 热传导方程差分格式的收敛性和稳定性在实际研究物理问题过程中, 往往能给出问题相应的数学表达式, 但是由于实际物理问题的复杂性, 它的解却一般不容易求出. 由此计算物理应运而生, 计算物理是以计算机为工具, 应用数学的方法解决物理问题的一门应用性学科, 是物理、数学和计算机三者结合的交叉性学科. 它产生于二战期间美国对核武器的研究, 伴随着计算机的发展而发展. 计算物理的目的不仅仅是计算, 而是要通过计算来解释和发现新的物理规律. 这一点它与传统的实验物理和理论物理并无差别, 所不同的只是使用的工具和方法. 计算物理早已与实验物理和理论物理形成三足鼎立之势, 甚至有人提出它将成为现代物理大厦的“栋梁”. 在一个物理问题中一个数值解往往比一个式子更直观, 更有价值. 在实际求解方程时, 除了一些特殊的情况下可以方便地求得其精确解外, 在一般情况下, 当方程或定解条件具有比较复杂的形式, 或求解区域具有比较复杂的形状时, 往往求不到, 或不易求到其精确解. 这就需要我们去寻找方程的近似解, 特别是数值近似解, 简称数值解. 这里主要研究的是热传导方程. 有限差分法是微分方程和积分微分方程数值解的方法. 其基本思想是把连续的定解区域用有限个离散点构成的网格来代替, 这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似, 于是原微分方程和定解条件就近似地代之以代数方程组, 即有限差分方程组, 解此方程组就可以得到原问题在离散点上的近似解. 然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解. 热传导的差分法是求解热传导方程的重要方法之一. 对于差分格式的的求解, 我们首先要关注差分格式的收敛性和稳定性. 对于一个微分方程建立的各种差分格式, 为了有实用意义, 一个基本要求是它们能够任意逼近微分方程, 即相容性要求. 一个差分格式是否有用, 就要看差分方程的精确解能否任意逼近微分方程的解, 即收敛性的概念. 此外, 还有一个重要的概念必须考虑, 即差分格式的稳定性. 因为差分格式的计

一维导热方程 有限差分法 matlab实现

第五次作业（前三题写在作业纸上） 一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf 文件，热扩散系数α=const ， 22T T t x α??=?? 1. 用Tylaor 展开法推导出FTCS 格式的差分方程 2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。 3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。 4. 编写M 文件求解上述方程，并用适当的文字对程序做出说明。（部分由网络搜索得到，添加，修改后得到。） function rechuandaopde %以下所用数据，除了t 的范围我根据题目要求取到了20000，其余均从pdf 中得来 a=0.00001;%a 的取值 xspan=[0 1];%x 的取值范围 tspan=[0 20000];%t 的取值范围 ngrid=[100 10];%分割的份数，前面的是t 轴的，后面的是x 轴的 f=@(x)0;%初值 g1=@(t)100;%边界条件一 g2=@(t)100;%边界条件二 [T,x,t]=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);%计算所调用的函数 [x,t]=meshgrid(x,t); mesh(x,t,T);%画图，并且把坐标轴名称改为x ，t ，T xlabel('x') ylabel('t') zlabel('T') T%输出温度矩阵 dt=tspan(2)/ngrid(1);%t 步长 h3000=3000/dt;

h9000=9000/dt; h15000=15000/dt;%3000,9000,15000下，温度分别在T矩阵的哪些行T3000=T(h3000,:) T9000=T(h9000,:) T15000=T(h15000,:)%输出三个时间下的温度分布 %不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面 %稳定性讨论,傅里叶级数法 dx=xspan(2)/ngrid(2);%x步长 sta=4*a*dt/(dx^2)*(sin(pi/2))^2; if sta>0,sta<2 fprintf('\n%s\n','有稳定性') else fprintf('\n%s\n','没有稳定性') error end %真实值计算 [xe,te,Te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid); [xe,te]=meshgrid(xe,te); mesh(xe,te,Te);%画图，并且把坐标轴名称改为xe，te，Te xlabel('xe') ylabel('te') zlabel('Te') Te%输出温度矩阵 %误差计算 jmax=1/dx+1;%网格点数 [rms]=wuchajisuan(T,Te,jmax) rms%输出误差

差分方程的解法分析及MATLAB实现(程序)

差分方程的解法分析及MATLAB 实现（程序） 摘自：张登奇,彭仕玉.差分方程的解法分析及其MATLAB 实现[J]. 湖南理工学院学报.2014(03) 引言 线性常系数差分方程是描述线性时不变离散时间系统的数学模型,求解差分方程是分析离散时间系统的重要内容.在《信号与系统》课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域 法[1]. 1 迭代法 例1 已知离散系统的差分方程为)1(3 1)()2(81)1(43)(-+=-+--n x n x n y n y n y ,激励信号为)()4 3()(n u n x n =,初始状态为21)2(4)1(=-=-y y ，.求系统响应. 根据激励信号和初始状态,手工依次迭代可算出24 59)1(,25)0(==y y . 利用MATLAB 中的filter 函数实现迭代过程的m 程序如下: clc;clear;format compact; a=[1,-3/4,1/8],b=[1,1/3,0], %输入差分方程系数向量,不足补0对齐 n=0:10;xn=(3/4).^n, %输入激励信号 zx=[0,0],zy=[4,12], %输入初始状态 zi=filtic(b,a,zy,zx),%计算等效初始条件 [yn,zf]=filter(b,a,xn,zi),%迭代计算输出和后段等效初始条件 2 时域经典法 用时域经典法求解差分方程:先求齐次解；再将激励信号代入方程右端化简得自由项,根据自由项形 式求特解；然后根据边界条件求完全解[3].用时域经典法求解例1的基本步骤如下. （1）求齐次解.特征方程为081432=+-αα,可算出4 1 , 2121==αα.高阶特征根可用MATLAB 的roots 函数计算.齐次解为. 0 , )4 1()21()(21≥+=n C C n y n n h （2）求方程的特解.将)()4 3()(n u n x n =代入差分方程右端得自由项为 ?????≥?==-?+-1,)4 3(9130 ,1)1()43(31)()43(1n n n u n u n n n 当1≥n 时,特解可设为n p D n y )4 3()(=,代入差分方程求得213=D . （3）利用边界条件求完全解.当n =0时迭代求出25)0(=y ,当n ≥1时,完全解的形式为 ,)4 3(213 )41()21()(21n n n C C n y ?++=选择求完全解系数的边界条件可参考文[4]选)1(),0(-y y .根据边界条件求得35,31721=-=C C .注意完全解的表达式只适于特解成立的n 取值范围,其他点要用 )(n δ及其延迟表示,如果其值符合表达式则可合并处理.差分方程的完全解为

有限差分法的Matlab程序(椭圆型方程)

有限差分法的Matlab程序（椭圆型方程） function FD_PDE(fun,gun,a,b,c,d) % 用有限差分法求解矩形域上的Poisson方程 tol=10^(-6); % 误差界 N=1000; % 最大迭代次数 n=20; % x轴方向的网格数 m=20; % y轴方向的网格数 h=(b-a)/n; % x轴方向的步长 l=(d-c)/m; % y轴方向的步长 for i=1:n-1 x(i)=a+i*h; end % 定义网格点坐标 for j=1:m-1 y(j)=c+j*l; end % 定义网格点坐标 u=zeros(n-1,m-1); %对u赋初值 % 下面定义几个参数 r=h^2/l^2; s=2*(1+r); k=1; % 应用Gauss-Seidel法求解差分方程 while k<=N % 对靠近上边界的网格点进行处理 % 对左上角的网格点进行处理 z=(-h^2*fun(x(1),y(m-1))+gun(a,y(m-1))+r*gun(x(1),d)+r*u(1,m-2)+u(2,m-1))/s; norm=abs(z-u(1,m-1)); u(1,m-1)=z; % 对靠近上边界的除第一点和最后点外网格点进行处理 for i=2:n-2 z=(-h^2*fun(x(i),y(m-1))+r*gun(x(i),d)+r*u(i,m-2)+u(i+1,m-1)+u(i-1,m-1))/s; if abs(u(i,m-1)-z)>norm; norm=abs(u(i,m-1)-z); end u(i,m-1)=z; end % 对右上角的网格点进行处理 z=(-h^2*fun(x(n-1),y(m-1))+gun(b,y(m-1))+r*gun(x(n-1),d)+r*u(n-1,m-2)+u(n-2,m-1))/s; if abs(u(n-1,m-1)-z)>norm norm=abs(u(n-1,m-1)-z); end u(n-1,m-1)=z; % 对不靠近上下边界的网格点进行处理 for j=m-2:-1:2 % 对靠近左边界的网格点进行处理

电磁场实验一_有限差分法的matlab实现

电磁场与电磁波实验报告 实验项目：_______有限差分法__ ____ 班级：_____ __12电子2 ____ __ 实验日期：__2014年12月23日 姓名：___ _ __陈奋裕 __ __ 学号：___ ___1215106003 _____ 组员姓名：___ _ __ __ __ 组员学号：___ ___ _____ 指导教师：_ ____张海 ______

一、实验目的及要求 1、学习有限差分法的原理与计算步骤； 2、学习用有限差分法解静电场中简单的二维静电场边值问题； 3、学习用Matlab 语言描述电磁场与电磁波中内容，用matlab 求解问题并用图形表示出了，学习matlab 语言在电磁波与电磁场中的编程思路。 二、实验内容 理论学习：学习静电场中边值问题的数值法中的优先差分法的求解知识； 实践学习：学习用matlab 语言编写有限差分法计算二维静电场边值问题； 三、实验仪器或软件 电脑（WIN7）、Matlab7.11 四、实验原理 基本思想是把连续的定解区域用有限个离散点构成的网格来代替， 这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似， 积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组 ， 解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 简单迭代法： 这一方法的求解过程是，先对场域内的节点赋予迭代初值(0),i j ?，这里上标(0)表示0次 （初始）近似值。然后按Laplace 方程 (k 1)(k)(k)(k)(k),1,,11,,11 []4 i j i j i j i j i j ?????+--++=+++(i,j=1,2,…) 进行反复迭代（k=0，1,2,…）。若当第N 次迭代以后，所有的内节点的相邻两次迭代值之间的最大误差不超过允许范围，即 (N)(N-1) ,,max|-|

用matlab实现线性常系数差分方程的求解

数字信号处理课程设计 题目：试实现线性常系数差分方程的求解 学院： 专业： 班级： 学号： 组员： 指导教师：

题目：用Matlab 实现线性常系数差分方程求解 一． 设计要求 1． 掌握线性常系数差分方程的求解 2． 熟练掌握Matlab 基本操作和各类函数调用 3． 结合Matlab 实现线性常系数差分方程的求解 二．设计原理 1．差分与差分方程 与连续时间信号的微分及积分运算相对应，离散时间信号有差分及序列求和运算。设有序列f(k),则称…,f(k+2),f(k+1),…,f(k －1),f(k －2),…为f(k)的移位序列。序列的差分可以分为前向差分和后向差分。一阶前向差分定义为 ()(1)()f k f k f k ?=+- （3.1—1） 一阶后向差分定义为 ()()(1)f k f k f k ?=-- （3.1—2） 式中Δ和Δ称为差分算子。由式（3.1—1）和式（3.1—2）可见，前向差分与后向差分的关系为 ()(1)f k f k ?=?- （3.1—3） 二者仅移位不同，没有原则上的差别，因而它们的性质也相同。此处主要采用后向差分，并简称其为差分。 由查分的定义，若有序列1()f k 、2()f k 和常数1a ，2a 则 1122112211221112221122[()()][()()][(1)(1)][()(1)][()(1)]()() a f k a f k a f k a f k a f k a f k a f k f k a f k f k a f k a f k ?+=+--+-=--+--=?+? （3.1—4） 这表明差分运算具有线性性质。 二阶差分可定义为 2()[()][()(1)]()(1) ()2(1)(2) f k f k f k f k f k f k f k f k f k ?=??=?--=?-?-=--+- （3.1—5） 类似的，可定义三阶、四阶、…、n 阶差分。一般地，n 阶差分

基于MATLAB的线性常系数差分方程求解

天津城市建设学院 课程设计任务书 2012 —2013 学年第 1 学期 计算机与信息工程 学院 电子信息工程 系 专业 课程设计名称： 数字信号处理 设计题目： 基于MATLAB 的线性常系数差分方程求解 完成期限：自2012 年 12月 17 日至 2012 年 12月 28 日共 2 周 设计依据、要求及主要内容： 一．课程设计依据 《数字信号处理》是电子信息类专业极其重要的一门专业基础课程，这门课程是将信号和系统抽象成离散的数学模型，并从数学分析的角度分别讨论信号、系统、信号经过系统、系统设计（主要是滤波器）等问题。 采用仿真可帮助学生加强理解，在掌握数字信号处理相关理论的基础上，根据数字信号处理课程所学知识，利用Matlab 对线性常系数方程进行求解，分析不同初始条件对解的影响。 二．课程设计内容 1、自行产生一个序列，对序列进行差分运算，并画出差分序列的时域波形图； 2、已知一个二阶线性常系数差分方程用下式表示 y(n)+a1y(n-1)+a2y(n-2)= b0x(n)+b1x(n-1)+b2x(n-2) 要求：（1）参数a1、a2、b0、b1、b2由运行时输入；（2）已知输入()0.5()n x n u n 画出x(n)的时域波形图；（3）求出x(n )的共轭对称分量e ()x n 和共轭反对称分量o ()x n ，并分别画出时域波形图；（4） 初始条件由运行时输入，求输出y(n)，并画出其波形；（5）对于不同的初始条件分析其输出是否一致，从中得出什么结论。 三．课程设计要求 1. 要求独立完成设计任务。 2. 课程设计说明书封面格式要求见《天津城市建设学院课程设计教学工作规范》附表1 3. 课程设计的说明书要求简洁、通顺，计算正确，图纸表达内容完整、清楚、规范。 4. 测试要求：根据题目的特点，编写Matlab 程序，绘制结果图形，并从理论上进行分析。 5. 课设说明书要求： 1) 说明题目的设计原理和思路、采用方法及设计流程。 2) 详细介绍运用的理论知识和主要的Matlab 程序。 3) 绘制结果图形并对仿真结果进行详细的分析。

用MATLAB仿真离散系统差分方程

HEFEI UNIVERSITY 信号与系统项目设计报告 系别电子信息与电子工程系 题目项目第十题 专业电子信息工程 班级 11电子信息工程（2）班 小组成员钟文俊（1105012012）谢伟明（1105012041)授课老师纪平 完成时间 2014.01.02

用MATLAB仿真离散系统差分方程 一、设计题目 -n f - + - n n n y y n = y f (- ( )2 ) ( )1 ( 5.0 .0 ) 25 ( )1 + 系统输入序列为) f nε =。 n (n )5.0( ) ( 二、设计要求 1、试用MATLAB绘出输入序列的时域波形； 2、用MATLAB求出该系统0～20区间的样值； 3、用MATLAB画出系统的零状态响应波形。 三、功能分析 差分方程反映的是关于离散变量的取值与变化规律。通过建立一个或几个离散变量取值所满足的平衡关系，从而建立差分方程。差分方程就是针对要解决的目标，引路系统或过程中的离散变量，根据实际背景的规律、性质、平衡关系，建立离散变量所满足的平衡关系等式，从而建立差分方程。通过求出和分析方程的解，或者分析得到方程解的特别性质（平衡性、稳定性、渐近性、振动性周期性等），从而把握这个离散变量的变化过程的规律，进步再结合其他分析，得到原问题的解。 四、设计原理分析 1、差分方程定义 含有未知函数yt=f(t)以及yt的差分Dyt， D2yt，…的函数方程，称为常差分方程(简称差分方程)；出现在差分方程中的差分的最高阶数，称为差分方程的阶。n阶差分方程的一般形式为F(t，yt，Dyt，…， Dnyt)=0，其中F是t，yt, Dyt，…， Dnyt的已知函数，且Dnyt一定要在方程中出现。 含有两个或两个以上函数值yt，yt+1，…的函数方程，称为(常)差分方程，出现在差分方程中未知函数下标的最大差，称为差分方程的阶。n阶差分方程的一般形式为F(t，yt，yt+1，…，yt+n)=0,其中F为t，yt，yt+1，…，yt+n的已知函数，且yt和yt+n一定要在差分方程中出现。 2、差分方程的意义与应用 差分方程模型有着广泛的应用。实际上，连续变量可以用离散变量来近似逼

怎么用matlab把传递函数转成差分方程

怎么用matlab把传递函数转成差分方程 以下是PID控制的部分代码（matlab的m文件）： ts=;采样时间= sys=tf(400,[1,50,0]);建立被控对象传递函数 dsys=c2d(sys,ts,'z');把传递函数离散化（问题1） [num,den]=tfdata(dsys,'v');离散化后提取分子、分母 rin=;输入为阶跃信号 u_1=; u_2=; 什么东西的初始状态（问题2） y_1=; y_2=; 是不是输出的初始状态 error_1=0;初始误差 x=[0 0 0]';PID的3个参数Kp Ki Kd组成的数组 p=100;仿真时间100ms for k=1:1:p r(k)=rin; u(k)=kpidi(1)*x(1)+kpidi(2)*x(2)+kpidi(3)*x(3) if u(k)>=10 u(k)=10; end if u(k)<=-10 u(k)=-10; end yout(k)=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2;（问题3） error(k)=r(k)-yout(k); %返回pid参数 u_2=u_1;u_1=u(k); y_2=y_1;y_1=yout(k); x(1)=error(k); x(2)=(error(k)-error_1)/ts; x(3)=x(3)+error(k)*ts; error_2=error_1; error_1=error(k); end 问题1：把传递函数离散化[SYSD,G]=C2D(SYSC,Ts,METHOD)这里面的method有好多种，而且用的method不一样得出的结果也不一样，这些参数究竟有什么区别（不要把matlab的help 给我翻译一遍，帮忙解释详细点）

用Matlab实现差分方程

（1）用Matlab实现差分方程 (i) y(n)=0.6y(n-1)- 0.08y(n-2)+x(n)； (ii) y(n)=0.7y(n-1)- 0.1y(n-2)+2x(n)-x(n-2)； 程序如下： B1=1; A1=[1,0.6,0.08]; xn=[1, zeros(1, 30)]; xi=filtic(B1,A1,xn); hn1=filter(B1,A1,xn,xi); n=0:length(hn1)-1; subplot(3,2,1);stem(n,hn1,'.') title('a'); xlabel('n');ylabel('h(n)') xn=ones(1,30); sn1=filter(B1,A1,xn,xi); n=0:length(sn1)-1; subplot(3,2,2);stem(n,sn1,'.') title('b'); xlabel('n');ylabel('s(n)') 系统2: B2=[2,0,-1]; A2=[1, -0.7,0.1]; xn=[1, zeros(1, 30)]; xi=filtic(B2,A2,xn); hn2=filter(B2,A2,xn,xi); n=0:length(hn1)-1; subplot(3,2,3);stem(n,hn2,'.') title('c'); xlabel('n');ylabel('h(n)') xn=ones(1,30); sn2=filter(B1,A1,xn,xi); n=0:length(sn1)-1; subplot(3,2,2);stem(n,sn2,'.') title('d'); xlabel('n');ylabel('s(n)') （2）设3项滑动平均滤波器的输入信号为sin(n/6)u(n) ,画出该滤波器的输入和输出的前15个序列值。 程序如下： B=[1,1,1]/3; A=1; n=0:30; xn=sin(n*pi/6); yn=filter(B,A,xn); subplot(3,2,1);stem(n, xn,'.') ; ylabel('x(n)'); subplot(3,2,3);stem(n, yn,'.') ; xlabel('n');ylabel('y(n)');

基于MATLAB的线性常系数差分方程求解.

数字信号处理课程设计 题目：基于MATLAB的线性常系数差分方程求解学院： 专业： 班级： 学号： 姓名： 指导教师：

目录 摘要 (1) 第一章背景 (3) 1.1 背景知识 (3) 1.2 《数字信号课程》特点 (3) 1.3 软件介绍 (4) 1.4 MATLAB及数字信号处理 (4) 第二章设计目的及要求 (6) 2.1 设计目的 (6) 2.2 课程设计的内容要求 (7) 2.2.1 设计要求 (7) 第三章设计任务 (8) 第四章设计原理 (9) 4.1 差分与差分方程 (9) 4.2 线性常系数差分方程 (14) 4.3 线性常系数差分方程的求解 (15) 第五章设计过程 (16) 5.1 用MATLAB求解差分方程 (16) 第六章设计代码及结果 (18) 6.1 MATLAB源程序 (18) 6.2 程序运行结果 (20) 6.3 比较结果总结 (24) 第七章收获与体会 (25) 致谢 (27) 参考文献 (28)

摘要 《数字信号处理》分析了数字信号处理课程的重要性及特点，为了帮助学生理解与掌握课程中的基本概念、基本原理、基本分析方法，提出了用MATLAB进行数字信号处理课程设计的思路，并阐述了课程设计的具体方法、步骤和内容。 MATLAB语言是一种广泛应用于工程计算及数值分析领域的新型高级语言，MATLAB功能强大、简单易学、编成效率高，深受广大科技工作者的喜爱，特别是MATLAB还具有信号分析工具箱，不需具备很强的编程能力，就可以很方便地进行语音信号分析、处理和设计。 线性常系数差分方程求解是数字信号处理课程中常出现的课题，也是现代科学中值得深入研究的一个课题 本文介绍了线性常系数差分方程的基本概念，论述了其求解方法，并用MATLAB具体实现了线性常系数差分方程的求解。 基于MATLAB的线性常系数差分方程求解主要是用MATLAB作为工具平台，设计中涉及到差分方程的递推求解以及用filter对系数向量的归一化等等。通过数字信号处理课程的理论知识的综合运用，从实践上初步实现对数字信号的处理。 关键字：MATLAB，线性常系数差分方程，数字信号处理。

matlab实现有限差分法计算电场强度(最新)

实验一：有限差分法研究静电场边值问题 实验报告人：年级和班级：学号： 1. 实验用软件工具： Matlab 2. 实验原理：电磁场课本P36-38 1)差分方程 2)差分方程组的解 简单迭代法 高斯-赛德尔迭代法 逐次超松弛法 3. 实验步骤： 1）简单迭代法 程序： hx=41;hy=21; v1=zeros(hy,hx); v1(hy,:)=zeros(1,hx); v1(1,:)=ones(1,hx)*100; v1(:,1)=zeros(hy,1); v1(:,hx)=zeros(hy,1); v1 v2=v1;maxt=1;t=0; k=0; while(maxt>1e-5) k=k+1; maxt=0; for i=2:hy-1 for j=2:hx-1 v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4; t=abs(v2(i,j)-v1(i,j)); if(t>maxt) maxt=t;end end end v1=v2; end v2 k clf subplot(1,2,1),mesh(v2) axis([0,41,0,21,0,100]) subplot(1,2,2),contour(v2,15) hold on

axis([-1,42,-1,25]) plot([1,1,hx,hx,1],[1,hy+1,hy+1,1,1],'r') text(hx/2,0.3,'0V','fontsize',11); text(hx/2-0.5,hy+0.5,'100V','fontsize',11); text(-0.5,hy/2,'0V','fontsize',11); text(hx+0.3,hy/2,'0V','fontsize',11); hold off 当W=1e-5, 迭代次数：1401次 2）高斯-赛德尔迭代法 程序： hx=41;hy=21; v1=ones(hy,hx); v1(hy,:)=zeros(1,hx); v1(1,:)=ones(1,hx)*100; v1(:,1)=zeros(hy,1); v1(:,hx)=zeros(hy,1); v2=v1;maxt=1;t=0; k=0; while(maxt>1e-5) k=k+1; maxt=0; for i=2:hy-1 for j=2:hx-1 v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4; t=abs(v2(i,j)-v1(i,j)); if(t>maxt) maxt=t;end end end v1=v2; end v2 k clf subplot(1,2,1),mesh(v2) axis([0,41,0,21,0,100]) subplot(1,2,2),contour(v2,15) hold on axis([-1,42,-1,25]) plot([1,1,hx,hx,1],[1,hy+1,hy+1,1,1],'r') text(hx/2,0.3,'0V','fontsize',11); text(hx/2-0.5,hy+0.5,'100V','fontsize',11); text(-0.5,hy/2,'0V','fontsize',11); text(hx+0.3,hy/2,'0V','fontsize',11); hold off

有限差分法解偏微分方程

有限差分法解偏微分方程综述 绪论 有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor 级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。 从差分的空间形式来考虑，可分为中心格式和逆风格式。 考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。 目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式， 目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。 有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。 有限差分法求解偏微分方程 在采用数值计算方法求解偏微分方程时，若将每一处导数由有限差分近似公式替代，从而把求解偏微分方程的问题转换成求解代数方程的问题，即所谓的有限差分法。有限差分法求解偏微分方程的步骤如下： 1、区域离散化，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格； 2、近似替代，即采用有限差分公式替代每一个格点的导数； 3、逼近求解。换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程 有限差分法的应用 抛物型方程的差分方法 1. 简单差分法

Matlab求解微分方程(组)及偏微分方程(组)

创作编号：BG7531400019813488897SX 创作者：别如克* 第四讲Matlab求解微分方程（组） 理论介绍：Matlab求解微分方程（组）命令 求解实例：Matlab求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介 1.在Matlab中，用大写字母D表示导数，Dy表示y关于自变量的一阶导数，D2y表示y关于自变量的二阶导数，依此类推.函数dsolve用来解决常微分方程（组）的求解问题，调用格式为： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解. 注意，系统缺省的自变量为t 2.函数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 =上从0 t 到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阶公式作误差估计来调节步长，具有中等的精度.

2.差分方程及其求解---数字信号处理实验报告

计算机与信息工程学院验证性实验报告 一、实验目的 1．学习并掌握系统的差分方程表示方法以及差分方程的相关概念。 2．熟练使用filter 函数对差分方程进行数值求解。 3．掌握差分方程的求解及MATLAB 实现方法。 二、实验原理及方法 1．一LTI 系统可以用一个线性常系数差分方程表示： 00()()N M k m k m a y n k b x n m ==-=-∑∑，任意n 如果N a ≠0，那么这个差分方程就是N 阶的，已知系统的输入序列，用这个方程可以根据当前输入x(n)和以前M 点的输入x(n-m ),…,x(n-1),以及以前N 点的输出y(n-N),…,y(n -1)来计算当前输出y(n)。在实际中这个方程在时间上是从n =-∞到n =+∞朝前计算的，因此该方程的另一种形式是： 00()()()M N m k m k y n b x n m a y n k ===---∑∑ 方程的解能以下面形式求得：()()()H p y n y n y n =+分别为方程的齐次解跟特解部分。已知输入和差分方程的稀疏，可用filter 对差分方程进行数值求解。最简单形式为：y=filter(b,a,x) 其中b=[b0,b1,…,bM];a=[a0,a1,…,aN]; 2．上面差分方程解的形式为齐次解和特解，另外还可以求零输入解和零状态解理论计算中要用到z 变换，请好好掌握z 变换的内容。用MATLAB 实现时，若已知初始条件，则应用y=filter(b,a, x, xic)来求完全响应。这里xic 是初始状态输入数组。MATLAB 还提供一种filtic 函数来得到xic 。

一维扩散方程的有限差分法matlab

一维扩散方程的有限差分法 ——计算物理实验作业七 陈万 ● 题目： 编程求解一维扩散方程的解 取1.0,1.0,1.0,10,0.1,0,1,1,0,1,1max 0222111======-=====τh D t a c b a c b a 。输出t=1,2,...,10时刻的x 和u(x),并与解析解u=exp(x+0.1t)作比较。 ● 主程序： % 一维扩散方程的有限差分法 clear,clc; %定义初始常量 a1 = 1; b1 = 1; c1 = 0; a2 = 1;b2 = -1; c2 = 0; a0 = 1.0; t_max = 10; D = 0.1; h = 0.1; tao = 0.1; %调用扩散方程子函数求解 u = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2); ● 子程序1： function output = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2) % 一维扩散方程的有限差分法，采用隐式六点差分格式(Crank-Nicolson) % a0: x 的最大值 % t:_max: t 的最大值 % h: 空间步长 % tao: 时间步长

% D：扩散系数 % a1,b1,c1是（x=0）边界条件的系数；a2,b2,c2是（x=a0）边界条件的系数x = 0:h:a0; n = length(x); t = 0:tao:t_max; k = length(t); P = tao * D/h^2; P1 = 1/P + 1; P2 = 1/P - 1; u = zeros(k,n); %初始条件 u(1,:) = exp(x); %求A矩阵的对角元素d d = zeros(1,n); d(1,1) = b1*P1+h*a1; d(2:(n-1),1) = 2*P1; d(n,1) = b2*P1+h*a2; %求A矩阵的对角元素下面一行元素e e = -ones(1,n-1); e(1,n-1) = -b2; %求A矩阵的对角元素上面一行元素f f = -ones(1,n-1); f(1,1) = -b1; R = zeros(k,n);%求R %追赶法求解