有限差分法求解偏微分方程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) 求解区域的网格划分步长参数如下：
第四讲 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 =上从0t 到f t
差分方程的解法分析及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 求解微分方程（组）命令 求解实例：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 =上从0t 到f t 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,, ,f t t t t 上的解，则令 tspan 012[,,,]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，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 阶差分
天津城市建设学院 课程设计任务书 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) 绘制结果图形并对仿真结果进行详细的分析。
（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)');
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把传递函数转成差分方程 以下是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) 第一章背景 (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，线性常系数差分方程，数字信号处理。
第五次作业（前三题写在作业纸上） 一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.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%输出误差
创作编号：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阶公式作误差估计来调节步长，具有中等的精度.
计算机与信息工程学院验证性实验报告 一、实验目的 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 。
南京理工大学 课程考核论文 课程名称：高等数值分析 论文题目：有限差分法求解偏微分方程姓名：罗晨 学号： 1 成绩： 有限差分法求解偏微分方程
一、主要内容 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 ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。 推导差分方程的过程中需要用到的泰勒展开公式如下： ()2100000000()()()()()()()......()(()) 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) 求解区域的网格划分步长参数如下： 11k k k k t t x x h τ ++-=?? -=? (2-2) 2.1 古典显格式
第十六章 偏微分方程的数值解法 科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。包含多个自变量的微分方程称为偏微分方程(partial differential equation)，简称PDE 。偏微分方程问题，其求解是十分困难的。除少数特殊情况外，绝大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。 16.1 几类偏微分方程的定解问题 一个偏微分方程的表示通常如下： (,,,,)xx xy yy x y A B C f x y Φ+Φ+Φ=ΦΦΦ (16.1.1) 式中，,,A B C 是常数，称为拟线性(quasilinear)数。通常，存在3种拟线性方程： 双曲型(hyperbolic)方程：240B AC ->； 抛物线型(parabolic)方程：240B AC -=； 椭圆型(ellliptic)方程：240B AC -<。 16.1.2 双曲型方程 最简单形式为一阶双曲型方程： 0u u a t x ??+=?? (16.1.2) 物理中常见的一维振动与波动问题可用二阶波动方程： 22222u u a t x ??=?? (16.1.3) 描述，它是双曲型方程的典型形式。方程的初值问题为：
222220 0,(,0)()()t u u a t x t x u x x u x x t ?ψ=???=>-∞<<+∞ ????? =?? ??=-∞<<+∞ ??? (16.1.4) 边界条件一般有三类，最简单的初边值问题为： 22222120 00,0(,0)(0,)(),(,)()0()t u u a t T x l t x u x l u t g t u l t g t t T u x x t ?ψ=???==<<<?? (16.1.8) 方程可以有两种不同类型的定解问题： (1) 初值问题： 2200,(,0)()u u a t x t x u x x x ????-=>-∞<<+∞? ????=-∞<<+∞ ? (16.1.6) (2) 初边值问题： 2212 00,0(,0)()0(0,)(),(,)()0u u a t T x l t x u x x x l u t g t u l t g t t T ????-=<<<
[原创]偏微分方程数值解法的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==7 C=1; end %计算步长 dx=uX/M;%x的步长 dt=uT/N;%t的步长 x=(0:M)*dx; t=(0:N)*dt;
实验二： 微分方程与差分方程模型Matlab 求解 一、实验目的 [1] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； [2] 熟悉MATLAB 软件关于微分方程求解的各种命令； [3] 通过范例学习建立微分方程方面的数学模型以及求解全过程； [4] 熟悉离散 Logistic 模型的求解与混沌的产生过程。 二、实验原理 1. 微分方程模型与MATLAB 求解 解析解 用MATLAB 命令dsolve(‘eqn1’,’eqn2’, ...) 求常微分方程（组）的解析解。其中‘eqni'表示第i 个微分方程，Dny 表示y 的n 阶导数,默认的自变量为t 。 （1） 微分方程 例1 求解一阶微分方程 21y dx dy += (1) 求通解 输入： dsolve('Dy=1+y^2') 输出： ans = tan(t+C1) （2）求特解 输入： dsolve('Dy=1+y^2','y(0)=1','x') 指定初值为1，自变量为x 输出： ans = tan(x+1/4*pi)
例2 求解二阶微分方程 221 ()0 4 (/2)2(/2)2/x y xy x y y y πππ '''++-=='=- 原方程两边都除以2x ，得211 (1)04y y y x x '''++-= 输入: dsolve('D2y+(1/x)*Dy+(1-1/4/x^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x') ans = - (exp(x*i)*(pi/2)^(1/2)*i)/x^(1/2) + (exp(x*i)*exp(-x*2*i)*(pi/2)^(3/2)*2*i)/(pi*x^(1/2)) 试试能不用用simplify 函数化简 输入: simplify(ans) ans = 2^(1/2)*pi^(1/2)/x^(1/2)*sin(x) （2）微分方程组 例3 求解 d f /d x =3f +4g ; d g /d x =-4f +3g 。 （1）通解: [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g') f = exp(3*t)*(C1*sin(4*t)+C2*cos(4*t)) g = exp(3*t)*(C1*cos(4*t)-C2*sin(4*t)) 特解： [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1') f = exp(3*t)*sin(4*t) g = exp(3*t)*cos(4*t)