油藏数值模拟显式差分MATLAB源程序

matlab实现数值计算功能源程序(个人整理)

matlab数值计算功能1,基础运算（1）多项式的创建与表达将多项式(x-6)(x-3)(x-8)表示为系数形式a=[6 3 8] % 写成根矢量pa=poly(a)% 求出系数矢量ppa=poly2sym(pa,'x') % 表示成符号形式ezplot(ppa,[-50,50])求3介方阵A的特征多项式a=[6 2 4;7 5 6;1 3 6 ];pa=poly(a)% 写出系数矢量ppa=poly2sym(pa) %表示成符号形式ezplot(ppa,[-50,50]) % 绘图求x^3-6x^2-72x-27的根。

a=[1,-6,-72,-85]; % 写出多项式系数矢量r=roots(a) % 求多项式的根（2）多项式的乘除运算c=conv(a,b) %乘法[q,r]=deconv(c,a)% 除法求a(s)=s^2+2s+3乘以b(s)=4s^2+5s+6的乘积a=[1 2 3]b=[4 5 6] % 写出系数矢量c=conv(a,b)c=poly2sym(c,'s') % 写成符号形式的多项式展开(s^2+2s+2)(s+4)(s+1)并验证结果被(s+4),(s+3)除后的结果。

c=conv([1,2,2],conv([1,4],[1,1]))cs=poly2sym(c,'s')c=[1 7 16 18 8][q1,r1]=deconv(c,[1,4])[q2,r2]=deconv(c,[1,3])cc=conv(q2,[1,3])test=((c-r2)==cc)其他常用的多项式运算命令pa=polyval(p,s) % 按数组规则计算给定s时多项式的值pm=polyvalm(p,s)% 按矩阵规则计算给定s时多项式的值[r,p,k]=residue(b,a) % 部分分式展开，b,a分别是分子，分母多项式系数矢量。

r,p,k分别是留数，极点和值项矢量。

matlab有限差分法

matlab有限差分法一、前言Matlab是一种广泛应用于科学计算和工程领域的计算机软件，它具有简单易学、功能强大、易于编程等优点。

有限差分法（Finite Difference Method）是一种常用的数值解法，它将微分方程转化为差分方程，通过对差分方程进行离散化求解，得到微分方程的数值解。

本文将介绍如何使用Matlab实现有限差分法。

二、有限差分法基础1. 有限差分法原理有限差分法是一种通过将微分方程转化为离散形式来求解微分方程的数值方法。

其基本思想是将求解区域进行网格划分，然后在每个网格点上进行逼近。

假设要求解一个二阶常微分方程：$$y''(x)=f(x,y(x),y'(x))$$则可以将其转化为离散形式：$$\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}=f(x_i,y_i,y'_i)$$其中$h$为网格步长，$y_i$表示在$x_i$处的函数值。

2. 一维情况下的有限差分法对于一维情况下的常微分方程：$$\frac{d^2 y}{dx^2}=f(x,y,y')$$可以使用中心差分法进行离散化：$$\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}=f(x_i,y_i,y'_i)$$这个方程可以写成矩阵形式：$$A\vec{y}=\vec{b}$$其中$A$为系数矩阵，$\vec{y}$为函数值向量，$\vec{b}$为右端项向量。

三、Matlab实现有限差分法1. 一维情况下的有限差分法假设要求解的方程为：$$\frac{d^2 y}{dx^2}=-\sin(x)$$首先需要确定求解区域和网格步长。

在本例中，我们将求解区域设为$[0,2\pi]$，网格步长$h=0.01$。

则可以通过以下代码生成网格：```matlabx = 0:0.01:2*pi;```接下来需要构造系数矩阵和右端项向量。

根据上面的公式，系数矩阵应该是一个三对角矩阵，可以通过以下代码生成：```matlabn = length(x)-2;A = spdiags([-ones(n,1), 2*ones(n,1), -ones(n,1)], [-1 0 1], n, n); ```其中`spdiags`函数用于生成一个稀疏矩阵。

有限差分法matlab程序

有限差分法matlab程序
有限差分法是一种常见的数值分析方法，它可以用来求解微分方程的解，称为数值解。

有限差分法的实现可以使用MATLAB，MATLAB 的强大的编程能力可以帮助我们更方便的实现有限差分法。

首先，我们需要了解什么是有限差分法。

有限差分法是一种数值计算的方法，它用一组等间距的数组来拟合一个函数，然后利用这些数组求出函数的值。

有限差分法可以用来解决复杂的微分方程，这是由于它能够将一个复杂的微分方程分解成一系列简单的微分方程，每个方程都可以用有限差分法来解决。

其次，使用MATLAB来实现有限差分法需要了解MATLAB中有限差分法的必要组件。

MATLAB中有限差分法的组件主要有：矩阵操作、函数计算、微分方程求解等。

首先，我们使用矩阵操作来实现函数的拟合，然后使用函数计算得到函数的值。

最后，我们使用微分方程求解器来解决复杂的微分方程。

最后，实现有限差分法在MATLAB中需要编写MATLAB程序。

MATLAB 程序首先需要编写矩阵操作程序，使用矩阵操作程序拟合函数；然后编写函数计算程序，计算函数的值；最后，编写微分方程求解程序，利用微分方程求解器来解决复杂的微分方程。

有限差分法是一种重要的数值分析方法，在很多科学中都有着广泛的应用，MATLAB可以让我们更方便的实现有限差分法，在微分方程的求解中发挥其独特的作用。

本文简要介绍了如何使用MATLAB来实现有限差分法，希望能为大家提供一些帮助。

matlab差分函数

matlab差分函数
Matlab差分函数是用于计算离散数据序列的差分。

差分是指通过对序列中的相邻元素进行减法运算，从而得到一个新序列的过程。

在Matlab中，使用diff函数可以快速计算序列的差分。

使用diff函数的基本语法如下：
diff(y)
其中，y为需要计算差分的序列。

diff函数返回一个新的序列，其元素值为相邻元素的差值。

例如，对于一个序列y=[1 3 5 7 9]，使用diff函数可以得到如下结果：
>> y = [1 3 5 7 9];
>> diff(y)
ans =
2 2 2 2
可以看到，diff函数返回了一个新的序列[2 2 2 2]，其元素值为相邻元素的差值。

除了基本语法外，diff函数还有一些常用的参数，例如：
diff(y,n)
其中，n表示需要进行差分的级数。

例如，如果n=2，则diff 函数将对y序列进行两次差分。

此外，diff函数还可以指定差分的方向，例如：
diff(y,[],1)
表示对y序列进行列方向的差分。

总之，Matlab差分函数可以快速计算序列的差分，是数据分析和处理的重要工具之一。

《油藏数值模拟》差分方程

二、差商
2、二阶差商
Pi+1
=
Pi
+
ΔxP' (
xi
)+
( Δx )2 2!
P'' (
xi
)+
( Δx )3 3!
P''' (
xi
)+
( Δx )4 4!
P( 4 )(
xi
)+⋅⋅⋅
Pi−1
=
Pi
− ΔxP' (
xi
)+
( Δx )2 2!
P'' (
xi
) − ( Δx )3 3!
P''' (
xi
)+
中国石油大学（北京）油藏数值模拟研究中心
第2节有限差分法
一、网格系统
1、全局正交网格(Globally Orthogonal Grid)
全局正交网格单元的外边界通常和坐标轴平行，而在外边界处，为了顺应边界的复杂形状，即可以对边界形状进行微弱的扭曲，也可以用台阶来近似。
一维网格
一维径向网格
中国石油大学（北京）油藏数值模拟研究中心
i-1
i
i+1
x0
xi-Δx
xi
xi+Δx
xN+1 x
中国石油大学（北京）油藏数值模拟研究中心
第2节有限差分法
二、差商
1、一阶差商
∂P
如图, 将P(x)在x点微商 ∂x xi ,表示为离散
Pi-1
Pi
Pi+1
x
i-1 i i+1
点上P(i=1,…,n)的线性函数。首先把 Pi+1 = P ( xi+1 ) 展开成x点泰勒级数:

油藏数值模拟隐式差分MATLAB源程序

油藏数值模拟隐式差分MATLAB源程序%隐式求解方法%t为投产后某一时刻，单位：天%d:迭代时间；%Pwf1:W1井底流压;%Q2:W2井产油量;function[P,d,Pwf1,Q2]=implict(t)%油藏参数Pini=20;u=5e-3;C=2e-4;Q1=30;Pwf2=15;d某=200;dy=200;dt=24;n=t某24/dt;%迭代时间步数re=0.208某d某;rw=0.1;%渗透率K=[02592222001901801850000;259259222200190180185185000;31031 0240235228210195195000;330330290270250230205197.51801850;3503503 00280259222200190180185185;340340320290310240235228210195195;355 355335315310290270250230205205;0000325300280240210215215;0000340 320290260235225225;00003553353152952752550];%厚度H=K/50;%孔隙度Fai=(K.某0.02+15)/100;%地层压力P=Pini某one(10,11);%P(1,:)=0;%P(2:10,1)=0;%P(2:10,11)=0;%P(2:4,8:10)=0;% P(8:10,2:4)=0;%系数矩阵%初始化a=zero(10,11);b=zero(10,11);c=zero(10,11);d=zero(10,11);e=ze ro(10,11);%fori=2:9forj=2:10a(i,j)=3600e-9某dt某2某H(i,j-1)某K(i,j-1)某K(i,j)/u/C/Fai(i,j)/d某/d某/(H(i,j-1)某K(i,j-1)+H(i,j)某K(i,j));b(i,j)=3600e-9某dt某2某H(i,j+1)某K(i,j+1)某K(i,j)/u/C/Fai(i,j)/d某/d某/(H(i,j+1)某K(i,j+1)+H(i,j)某K(i,j));c(i,j)=3600e-9某dt某2某H(i+1,j)某K(i+1,j)某K(i,j)/u/C/Fai(i,j)/d某/d某/(H(i+1,j)某K(i+1,j)+H(i,j)某K(i,j));d(i,j)=3600e-9某dt某2某H(i-1,j)某K(i-1,j)某K(i,j)/u/C/Fai(i,j)/d某/d某/(H(i-1,j)某K(i-1,j)+H(i,j)某K(i,j));e(i,j)=-1-a(i,j)-b(i,j)-c(i,j)-d(i,j);endend%初始时刻地层压力p=Pini某one(52,1);%系数矩阵AA=zero(52);%第9行A(1,1)=e(9,6);A(1,2)=b(9,6);A(1,6)=d(9,6);fori=2:4A(i,i-1)=a(9,i+5);A(i,i)=e(9,i+5);A(i,i+1)=b(9,i+5);A(i,i+5)=d(9,i+5); endA(5,4)=a(9,10);A(5,5)=e(9,10)+b(9,10);A(5,10)=d(9,10);%第8行A(6,1)=c(8,6);A(6,6)=e(8,6);A(6,7)=b(8,6);A(6,12)=d(8,6);fori=7:9A(i,i-5)=c(8,i);A(i,i-1)=a(8,i);A(i,i)=e(8,i);A(i,i+1)=b(8,i);A(i,i+6)=d(8,i);endA(10,5)=c(8,10);A(10,9)=a(8,10);A(10,10)=e(8,10)+b(8,10);A(10,16)=d(8,10);%第7行A(11,11)=e(7,5);A(11,12)=b(7,5);A(11,20)=d(7,5);fori=12:15A(i,i-6)=c(7,i-6);A(i,i-1)=a(7,i-6);A(i,i)=e(7,i-6);A(i,i+1)=b(7,i-6);A(i,i+9)=d(7,i-6);endA(16,10)=c(7,10);A(16,15)=a(7,10);A(16,16)=e(7,10)+b(7,10);A(16,25)=d(7,10);%第6行A(17,17)=e(6,2)+a(6,2);A(17,18)=b(6,2);A(17,26)=d(6,2); fori=18:19A(i,i-1)=a(6,i-15);A(i,i)=e(6,i-15);A(i,i+1)=b(6,i-15);A(i,i+9)=d(6,i-15);endfori=20:23A(i,i-9)=c(6,i-15);A(i,i-1)=a(6,i-15);A(i,i)=e(6,i-15);A(i,i+1)=b(6,i-15);A(i,i+9)=d(6,i-15);end%w2井定井底流压生产A(24,15)=c(6,9);A(24,23)=a(6,9);A(24,24)=e(6,9)-(3600e-9)某2某pi某K(6,9)某dt/u/C/Fai(6,9)/d某/dy/log(re/rw);A(24,25)=b(6,9);A(24,33)=d(6,9);A(25,16)=c(6,10);A(25,24)=a(6,10);A(25,25)=e(6,10)+b(6,10);A(25,34)=d(6,10);%第5行A(26,17)=c(5,2);A(26,26)=e(5,2)+a(5,2);A(26,27)=b(5,2);A(26,35)=d(5,2); fori=27:31A(i,i-9)=c(5,i-24);A(i,i-1)=a(5,i-24);A(i,i)=e(5,i-24);A(i,i+1)=b(5,i-24);A(i,i+9)=d(5,i-24);endfori=32:33A(i,i-9)=c(5,i-24);A(i,i-1)=a(5,i-24);A(i,i)=e(5,i-24)+d(5,i-24);A(i,i+1)=b(5,i-24);endA(34,25)=c(5,10);A(34,33)=a(5,10);A(34,34)=e(5,10)+b(5,10)+d(5,10);%第4行A(35,26)=c(4,2);A(35,35)=e(4,2)+a(4,2);A(35,36)=b(4,2);A(35,41)=d(4,2);fori=36:39A(i,i-9)=c(4,i-33);A(i,i-1)=a(4,i-33);A(i,i)=e(4,i-33);A(i,i+1)=b(4,i-33);A(i,i+6)=d(4,i-33);endA(40,31)=c(4,7);A(40,39)=a(4,7);A(40,40)=e(4,7)+b(4,7);A(40,46)=d(4,7);%第3行A(41,35)=c(3,2);A(41,41)=e(3,2)+a(3,2);A(41,42)=b(3,2);A(41,47)=d(3,2);fori=42:45A(i,i-6)=c(3,i-39);A(i,i-1)=a(3,i-39);A(i,i)=e(3,i-39);A(i,i+1)=b(3,i-39);A(i,i+6)=d(3,i-39);endA(46,40)=c(3,7);A(46,45)=a(3,7);A(46,46)=e(3,7)+b(3,7);A(46,52)=d(3,7);%第2行A(47,41)=c(2,2);A(47,47)=e(2,2)+a(2,2)+d(2,2);A(47,48)=b(2,2);fori=48:51A(i,i-6)=c(2,i-45);A(i,i-1)=a(2,i-45);A(i,i)=e(2,i-45)+d(2,i-45);A(i,i+1)=b(2,i-45);end。

matlab差分编码

matlab差分编码
在MATLAB中进行差分编码可以通过以下步骤完成：
1. 初始化：首先，你需要定义你的初始值，这将是你差分编码的起始点。

2. 计算差分：然后，你需要计算当前值与前一个值的差。

这可以通过简单的减法操作完成。

3. 编码：将差分值进行编码，可以使用不同的编码方式，例如使用霍夫曼编码，或者使用其他更简单的编码方式。

下面是一个简单的示例代码，展示了如何在MATLAB中进行差分编码：% 初始化
data = [1, 2, 3, 4, 5]; % 原始数据
diff_data = []; % 差分数据
% 计算差分
for i = 2:length(data)
diff_data(i) = data(i) - data(i-1);
end
% 显示差分数据
disp(diff_data);
以上代码首先定义了一个原始数据数组data，然后定义了一个空数组diff_data用于存储差分值。

接着，使用一个for循环来计算每一个数据点与前一
个数据点之间的差值，并将结果存储在diff_data中。

最后，使用disp函数显示差分数据。

请注意，这只是一个非常基础的示例。

在实际应用中，你可能需要使用更复杂的编码技术，并且可能需要处理各种边界条件和错误情况。

matlab 有限差分法代码

有限差分法是一种数值方法，用于求解偏微分方程。

下面是一个使用MATLAB 实现有限差分法的简单示例。

这个示例使用有限差分法求解一维热传导方程：
matlab
% 参数设定
L = 10; % 空间长度
T = 1; % 时间长度
dx = 0.1; % 空间步长
dt = 0.01; % 时间步长
Nx = round(L/dx); % 空间网格点数
Nt = round(T/dt); % 时间网格点数
% 初始化
u = zeros(1, Nx); % 温度数组
u_old = u; % 上一步的温度数组
% 初始条件
u(1) = 1; % 初始温度在 x=0 处为 1
u(Nx) = 0; % 初始温度在 x=L 处为 0
% 时间循环
for n = 1:Nt
% 空间循环
for i = 2:Nx-1
u_old(i) = u(i); % 保存当前温度值
u(i) = u_old(i) + dt/dx*(u_old(i+1) - u_old(i)); % 有限差分法更新温度值
end
% 时间循环继续，更新下一步的温度值
end
% 绘图
plot(0:dx:L, u);
xlabel('x');
ylabel('u');。

中心差分法 matlab代码

中心差分法(matlab代码)中心差分法是一种常用的数值求导方法，它利用函数在一点的两侧点进行逼近求导。

在matlab中，可以通过编写简单的代码来实现中心差分法的数值求导。

下面我将介绍如何使用matlab编写中心差分法的求导代码。

1. 准备工作在编写中心差分法的代码前，首先需要准备工作。

确保已经安装了matlab软件，并且已经打开了matlab编辑器。

需要确定要求导的函数，以及求导点的位置。

2. 编写函数在matlab中，可以使用函数来表示要求导的函数。

假设要求导的函数为f(x)，则可以使用如下代码来定义这个函数：```matlabfunction y = f(x)y = x^2; 示例：定义要求导的函数为x^2end```在这个示例中，我们定义了一个简单的函数f(x) = x^2作为要求导的函数。

3. 编写中心差分法代码编写中心差分法的代码需要考虑到求导点的选择。

中心差分法的原理是利用函数在求导点两侧的函数值来逼近求导值。

假设求导点为x0，假设步长为h，则中心差分法的求导公式为：```latexf'(x0) ≈ (f(x0+h) - f(x0-h)) / (2*h)```可以使用如下matlab代码来实现中心差分法的数值求导：```matlabfunction y = central_difference(x0, h)y = (f(x0 + h) - f(x0 - h)) / (2 * h);end```在这个示例中，我们定义了一个名为central_difference的函数，它接受两个参数x0和h，分别表示求导点的位置和步长。

在函数内部，我们使用了中心差分法的公式来计算数值导数的近似值。

4. 调用函数编写完中心差分法的代码后，可以通过调用这个函数来得到数值导数的近似值。

假设我们要在x=2的位置求函数f(x)=x^2的导数近似值，可以使用如下代码来进行计算：```matlabx0 = 2; 求导点的位置h = 0.01; 步长result = central_difference(x0, h); 调用central_difference函数进行计算disp(['数值导数的近似值为：', num2str(result)]); 显示计算结果```在这个示例中，我们通过调用central_difference函数来计算在x=2的位置的函数f(x)=x^2的导数近似值，并使用disp函数来显示计算结果。

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

差分方程的解法分析及MATLAB 实现（程序）摘自：张登奇,彭仕玉.差分方程的解法分析及其MATLAB 实现[J]. 湖南理工学院学报.2014(03) 引言线性常系数差分方程是描述线性时不变离散时间系统的数学模型,求解差分方程是分析离散时间系统的重要内容.在《信号与系统》课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域法[1].1 迭代法例1 已知离散系统的差分方程为)1(31)()2(81)1(43)(-+=-+--n x n x n y n y n y ,激励信号为)()43()(n u n x n =,初始状态为21)2(4)1(=-=-y y ，.求系统响应. 根据激励信号和初始状态,手工依次迭代可算出2459)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=+-αα,可算出41 , 2121==αα.高阶特征根可用MATLAB 的roots 函数计算.齐次解为. 0 , )41()21()(21≥+=n C C n y n n h （2）求方程的特解.将)()43()(n u n x n =代入差分方程右端得自由项为 ⎪⎩⎪⎨⎧≥⋅==-⋅+-1,)43(9130 ,1)1()43(31)()43(1n n n u n u n n n 当1≥n 时,特解可设为n p D n y )43()(=,代入差分方程求得213=D . （3）利用边界条件求完全解.当n =0时迭代求出25)0(=y ,当n ≥1时,完全解的形式为 ,)43(213 )41()21()(21n n n C C n y ⋅++=选择求完全解系数的边界条件可参考文[4]选)1(),0(-y y .根据边界条件求得35,31721=-=C C .注意完全解的表达式只适于特解成立的n 取值范围,其他点要用)(n δ及其延迟表示,如果其值符合表达式则可合并处理.差分方程的完全解为)(])43(213 )41(35)21(317[)1(])43(213 )41(35)21(317[)(25)(n u n u n n y n n n n n n ⋅+⋅+⋅-=-⋅+⋅+⋅-+=δ MATLAB 没有专用的差分方程求解函数,但可调用maple 符号运算工具箱中的rsolve 函数实现[5],格式为y=maple('rsolve({equs, inis},y(n))'),其中:equs 为差分方程表达式, inis 为边界条件,y(n)为差分方程中的输出函数式.rsolve 的其他格式可通过mhelp rsolve 命令了解.在MATLAB 中用时域经典法求解例1中的全响应和单位样值响应的程序如下.clc;clear;format compact;yn=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=(3/4)^n+1/3*(3/4)^(n-1),y(0)=5/2,y(-1)=4},y(n))'),hn=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=0,y(0)=1,y(1)=13/12},y(n))'),3 双零法根据双零响应的定义,按时域经典法的求解步骤可分别求出零输入响应和零状态响应.理解了双零法的求解原理和步骤,实际计算可调用rsolve 函数实现.yzi=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=0,y(-1)=4, y(-2)=12},y(n))'),yzs=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=(3/4)^n+1/3*(3/4)^(n-1),y(0)=1,y(-1)=0},y(n))'),4 变换域法设差分方程的一般形式为)()(00r n x b k n y a r Mr k N k -=-∑∑==.对差分方程两边取单边z 变换,并利用z 变换的位移公式得])()([])()([1010m r m r r M r l k l k k N k z m x z X z b z l y z Y z a ---=-=---=-=∑∑∑∑+=+整理成)()()()()()(00z X z X z B z Y z Y z A +=+形式有. )(, )(110110M M N N z b z b b z B z a z a a z A ----+++=+++=. )()(, )()(110110∑∑∑∑=--=--=--=--==M r r m m r r N k k l l k k z m x b s X zl y a s Y可以看出,由差分方程可直接写出 )(z A 和 )(z B ,系统函数)(/)()(z A z B z H =,将系统函数进行逆z 变换可得单位样值响应.由差分方程的初始状态可算出 )(0z Y ,由激励信号的初始状态可算出 )(0z X ,将激励信号进行z 变换可得 )(z X ,求解z 域代数方程可得输出信号的象函数 , )()()()()()(00z A z Y z X z X z B z Y -+= 对输出象函数进行逆z 变换可得输出信号的原函数)(n y .利用z 变换求解差分方程各响应的步骤可归纳如下:（1）根据差分方程直接写出 )(z A 、 )(z B 和)(z H ,)(z H 的逆变换即为单位样值响应；（2）根据激励信号算出 )(z X ,如激励不是因果序列则还要算出前M 个初始状态值；（3）根据差分方程的初始状态 )(, ),2( ),1(N y y y -⋅⋅⋅--和激励信号的初始状态 )(, ),2( ),1(M x x x -⋅⋅⋅--算出 )(0z Y 和 )(0z X ；（4）在z 域求解代数方程)()()()()()(00z X z X z B z Y z Y z A +=+得输出象函数 )(z Y , )(z Y 的逆变换即为全响应；（5）分析响应象函数的极点来源及在z 平面中的位置,确定自由响应与强迫响应,或瞬态响应与稳态响应；（6）根据零输入响应和零状态响应的定义,在z 域求解双零响应的象函数,对双零响应的象函数进行逆z 变换,得零输入响应和零状态响应.用变换域法求解例1的基本过程如下. 根据差分方程直接写出2181431 )(--+-=z z z A ,1311 )(-+=z z B .系统函数的极点为41,21. 对激励信号进行z 变换得)43/( )(-=z z z X .激励象函数的极点为3/4. 根据差分方程的初始状态算出102123 )(-+-=z z Y .根据激励信号的初始状态算出 0)(0=z X . 对z 域代数方程求解,得全响应的象函数)323161123/()83243125( )(2323-+-+-=z z z z z z z Y . 进行逆z 变换得全响应为)(])43(213 )41(35)21(317[)(n u n y n n n ⋅+⋅+⋅-= 其中,与系统函数的极点对应的是自由响应;与激励象函数的极点对应的是强迫响应. )(z Y 的极点都在z 平面的单位圆内故都是瞬态响应.零输入响应和零状态响应可按定义参照求解.上述求解过程可借助MATLAB 的符号运算编程实现.实现变换域法求解差分方程的m 程序如下: clc;clear;format compact;syms z n %定义符号对象% 输入差分方程、初始状态和激励信号%a=[1,-3/4,1/8],b=[1,1/3], %输入差分方程系数向量y0=[4,12],x0=[0], %输入初始状态,长度分别比a 、b 短1,长度为0时用[]xn=(3/4)^n, %输入激励信号,自动单边处理,u(n)可用1^n 表示% 下面是变换域法求解差分方程的通用程序，极点为有理数时有解析式输出 %N=length(a)-1;M=length(b)-1;%计算长度Az=poly2sym(a,'z')/z^N;Bz=poly2sym(b,'z')/z^M;%计算A(z)和B(z)Hz=Bz/Az;disp('系统函数H(z):'),sys=filt(b,a),%计算并显示系统函数hn=iztrans(Hz);disp('单位样值响应h(n)='),pretty(hn),%计算并显示单位样值响应Hzp=roots(a);disp('系统极点:');Hzp,%计算并显示系统极点Xz=ztrans(xn);disp('激励象函数X(z)='),pretty(Xz),%激励信号的单边z 变换Y0z=0;%初始化Y0(z),求Y0(z)注意系数标号与变量下标的关系for k=1:N;for l=-k:-1;Y0z = Y0z+a(k+1)*y0(-l)*z^(-k-l);endenddisp('初始Y0(z)'),Y0z,%系统初始状态的z 变换X0z=0;%初始化X0(z),求X0(z)注意系数标号与变量下标的关系for r=1:M;for m=-r:-1;X0z = X0z+b(r+1)*x0(-m)*z^(-r-m);endenddisp('初始X0(z)'),X0z,%激励信号起始状态的z 变换Yz=(Bz*Xz+X0z-Y0z)/Az;disp('全响应的z 变换Y(z)'),pretty(simple(Yz)),yn=iztrans(Yz);disp('全响应y(n)='),pretty(yn),% 计算并显示全响应Yziz=-Y0z/Az;disp('零输入象函数Yzi(z)='),pretty(Yziz),%零激励响应的z 变换yzin=iztrans(Yziz);disp('零输入响应yzi(n)='),pretty(yzin),% 计算并显示零输入响应 Yzsz=(Bz*Xz+X0z)/Az;disp('零状态象函数Yzs(z)='),pretty(Yzsz),%零状态响应的z 变换yzsn=iztrans(Yzsz);disp('零状态响应yzs(n)='),pretty(yzsn),% 计算并显示零状态响应该程序的运行过程与手算过程对应,显示在命令窗的运行结果与手算结果相同.。

matlab用差分法求解边值问题

MATLAB是一种强大的数学软件工具，用于解决各种数学问题，包括边值问题。

差分法是一种常用的数值分析方法，可以用于求解微分方程的边值问题。

本文将介绍如何使用MATLAB的差分法求解边值问题，包括算法原理、具体步骤和示例代码。

1. 算法原理差分法是一种数值解微分方程的常用方法，它的基本思想是用离散点代替连续函数，然后利用差分代替导数，将微分方程转化为代数方程，最终得到离散点上的数值解。

对于边值问题，差分法可以通过离散化区域，并在区域边界处使用相应的差分格式来求解。

2. 具体步骤（1）离散化区域需要将求解的区域进行离散化，将其划分为若干个离散点，形成一个网格。

这可以通过在区域内部建立均匀或非均匀的网格来实现。

（2）建立代数方程利用差分公式将原微分方程转化为代数方程。

根据边值问题的具体条件，可以构建相应的代数方程组。

（3）求解代数方程组利用MATLAB的线性方程求解器，如“backslash”运算符（\），求解得到代数方程组的解，即为边值问题的数值解。

3. 示例代码下面是一个简单的边值问题求解的MATLAB示例代码，以二阶常微分方程为例：``` MATLAB定义区域范围和离散步长a = 0;b = 1;N = 10; 离散点个数h = (b - a) / N; 离散步长构建代数方程组A = zeros(N-1, N-1);b = zeros(N-1, 1);for i = 1:N-1A(i, i) = 2/h^2;if i > 1A(i, i-1) = -1/h^2;endif i < N-1A(i, i+1) = -1/h^2;endb(i) = f(a + i*h); 右端项end求解代数方程组u = A \ b;绘制数值解x = a:h:b;plot(x, [0; u; 0]);```在这个示例代码中，我们首先定义了求解区域范围和离散步长，然后根据差分格式建立了代数方程组，并使用MATLAB的线性方程求解器求解得到数值解，最后绘制了数值解的图像。

matlab差分法

matlab差分法差分法（also known as finite difference method）是一种常见的数值计算方法，它将连续函数转化为离散的数值计算问题。

在MATLAB中，差分法可以用于求解微分方程、数值积分和数据插值等问题。

差分法的基本思想是通过近似表示函数的导数或积分，将连续函数转化为离散的数值计算问题。

差分的计算方式通常有前向差分、中心差分和后向差分三种。

前向差分是通过函数值在当前点和后一个点之间的差分来近似表示导数，计算公式如下：f'(x) ≈ (f(x+h) - f(x)) / h其中h表示插值点之间的步长，取值越小近似越精确，但计算量也会增加。

中心差分是通过函数值在前一个点和后一个点之间的差分来近似表示导数，计算公式如下：f'(x) ≈ (f(x+h) - f(x-h)) / (2h)中心差分法的精度比前向差分法高一阶，但计算量也相应增加。

后向差分是通过函数值在当前点和前一个点之间的差分来近似表示导数，计算公式如下：f'(x) ≈ (f(x) - f(x-h)) / h后向差分法的精度相对较低，但在一些特殊情况下仍然被广泛使用。

除了求导数，差分法还可以用于数值积分。

数值积分是通过对函数在一段区间上的离散点进行求和来近似表示整个区间上的积分值。

常见的数值积分方法包括矩形法、梯形法和辛普森法。

这些方法都可以通过对函数值使用相应的差分公式，来将积分转化为累加求和的形式，从而得到数值结果。

此外，差分法还可以用于数据插值。

数据插值是通过已知离散数据点的函数值，来计算未知点上的函数值。

常用的插值方法有线性插值、二次插值和三次插值等。

这些方法都可以通过使用差分公式，将插值问题转化为求解系数的线性方程组，进而得到插值结果。

综上所述，差分法是一种常见的数值计算方法，可以用于求解微分方程、数值积分和数据插值等问题。

MATLAB提供了丰富的函数和工具箱，使得差分法的实现变得简单和高效。

差分方程求解matlab代码

差分方程求解matlab代码差分方程在数值分析中是很重要的应用，matlab可以帮助我们有效的求解差分方程。

首先，我们要明确待求解的差分方程。

比如研究一个某种物质随时间的变化，可以用分类微分方程去描述：dC/dt = k * C其中C表示物质的浓度，t表示时间，k表示一个常量，dC/dt表示物质的瞬时浓度。

要使用matlab求解上面这个差分方程，需要如下步骤：1. 将差分方程写成matlab可以求解的形式。

因为 matlab 提供了 ode45 函数来求解常微分方程，所以我们需要将上面的差分方程转换成对应的常微分方程：dC/dt = k * C可以得到：dC/C = k * dt又即有：dC/C = k * dt2. 使用matlab编写代码进行求解。

假设已知t0=0, C0=3，t = 10（s），k = 2（1/s），则matlab代码如下：t0=0;t1=10; %定义起始时间、终止时间C0=3; %定义起始浓度k=2; %定义常量f=inline('k*x','x','k'); %定义右端函数[t,c]=ode45(f,[t0,t1],C0); %调用ode45函数进行求解plot(t,c), grid on; xlabel('Time(s)'); ylabel('Concentration(M)'); %绘制图像最后，使用 plot 函数绘制出解曲线，看到在 t0=0 时，物质的浓度为3 ，在 t=10 时，物质的浓度为 24。

在这里，我们用matlab求解了简单的一维差分方程，可以看到matlab还提供很多用于求解差分方程的函数，比如 ode23， ode113 等，这些函数可以帮助我们在不同的问题中有效地求解差分方程。

油藏数值模拟显式差分MATLAB源程序

%显式求解方法%t为投产后某一时刻，单位：天;%d:迭代时间；%Pwf1:W1井底流压;%Q2:W2井产油量;function [P1,d,Pwf1,Q2]=explict(t)%油藏参数Pini=20;u=5e-3;C=2e-4;Q1=30;Pwf2=15;dx=200;dy=200;dt=1.7;n=t*24/dt;%迭代时间步数re=0.208*dx;rw=0.1;%渗透率K=[0 259 222 200 190 180 185 0 0 0 0;259 259 222 200 190 180 185 185 0 0 0;310 310 240 235 228 210 195 195 0 0 0;330 330 290 270 250 230 205 197.5 180 185 0;350 350 300 280 259 222 200 190 180 185 185;340 340 320 290 310 240 235 228 210 195 195;355 355 335 315 310 290 270 250 230 205 205;0 0 0 0 325 300 280 240 210 215 215;0 0 0 0 340 320 290 260 235 225 225;0 0 0 0 355 335 315 295 275 255 0];%厚度H=K/50;%孔隙度Fai=(K.*0.02+15)/100;%原始地层压力P=Pini*ones(10,11);%P(8:10,1:4)=0;%P(1:3,9:11)=0;%P(1,1)=0;%P(1,8)=0;%P(10,11)=0;%P(4,11)=0;%生产后某一时刻地层压力P1=Pini*ones(10,11);%P1(8:10,1:4)=0;%P1(1:3,9:11)=0;%P1(1,1)=0;%P1(1,8)=0;%P1(10,11)=0;%P1(4,11)=0;%系数矩阵%初始化a=zeros(10,11);b=zeros(10,11);c=zeros(10,11);d=zeros(10,11);e=zeros(10,11);%for i=2:9for j=2:10a(i,j)=3600e-9*dt*2*H(i,j-1)*K(i,j-1)*K(i,j)/u/C/Fai(i,j)/dx/dx/(H(i,j-1)*K(i,j-1)+H(i,j)*K(i,j));b(i,j)=3600e-9*dt*2*H(i,j+1)*K(i,j+1)*K(i,j)/u/C/Fai(i,j)/dx/dx/(H(i,j+1)*K(i,j+1)+H(i,j)*K(i,j));c(i,j)=3600e-9*dt*2*H(i+1,j)*K(i+1,j)*K(i,j)/u/C/Fai(i,j)/dx/dx/(H(i+1,j)*K(i+1,j)+H(i,j)*K(i,j));d(i,j)=3600e-9*dt*2*H(i-1,j)*K(i-1,j)*K(i,j)/u/C/Fai(i,j)/dx/dx/(H(i-1,j)*K(i-1,j)+H(i,j)*K(i,j));e(i,j)=1-a(i,j)-b(i,j)-c(i,j)-d(i,j);endend% 循环次数tntn=1;while(tn<n)% 时间步长内压力变化量dPdP=zeros(10,11);%封闭边界压力计算for i=3:6P1(i,2)=c(i,2)*P(i+1,2)+(a(i,2)+e(i,2))*P(i,2)+b(i,2)*P(i,3)+d(i,2)*P(i-1,2);dP(i,2)= P1(i,2)-P(i,2);P1(i,1)=P1(i,2);endP1(2,2)=c(2,2)*P(3,2)+(a(2,2)+e(2,2)+d(2,2))*P(2,2)+b(2,2)*P(2,3);dP(2,2)= P1(2,2)-P(2,2);P1(2,1)=P1(2,2);P1(1,2)=P1(2,2);for j=3:6P1(2,j)=c(2,j)*P(3,j)+a(2,j)*P(2,j-1)+(e(2,j)+d(2,j))*P(2,j)+b(2,j)*P(2,j+1);dP(2,j)= P1(2,j)-P(2,j);P1(1,j)=P1(2,j);endP1(2,7)=c(2,7)*P(3,7)+a(2,7)*P(2,6)+(e(2,7)+b(2,7)+d(2,7))*P(2,7);dP(2,7)= P1(2,7)-P(2,7);P1(1,7)=P1(2,7);P1(2,8)=P1(2,7);for i=3:4P1(i,7)=c(i,7)*P(i+1,7)+a(i,7)*P(i,6)+(e(i,7)+b(i,7))*P(i,7)+d(i,7)*P(i-1,7);dP(i,7)= P1(i,7)-P(i,7);P1(i,8)= P1(i,7);endfor j=8:9P1(5,j)=c(5,j)*P(6,j)+a(5,j)*P(5,j-1)+(e(5,j)+d(5,j))*P(5,j)+b(5,j)*P(5,j+1);dP(5,j)= P1(5,j)-P(5,j);P1(4,j)= P1(5,j);endP1(5,10)=c(5,10)*P(6,10)+a(5,10)*P(5,9)+(e(5,10)+b(5,10)+d(5,10))*P(5,10);dP(5,10)= P1(5,10)-P(5,10);P1(4,10)=P1(5,10);P1(5,11)=P1(5,10);for i=6:9P1(i,10)=c(i,10)*P(i+1,10)+a(i,10)*P(i,9)+(e(i,10)+b(i,10))*P(i,10)+d(i,10)*P(i-1,10);dP(i,10)= P1(i,10)-P(i,10);P1(i,11)=P1(i,10);end%内部网格压力计算for i=8:9for j=6:9P1(i,j)=c(i,j)*P(i+1,j)+a(i,j)*P(i,j-1)+e(i,j)*P(i,j)+b(i,j)*P(i,j+1)+d(i,j)*P(i-1,j);dP(i,j)= P1(i,j)-P(i,j);endendfor j=5:9P1(7,j)=c(7,j)*P(8,j)+a(7,j)*P(7,j-1)+e(7,j)*P(7,j)+b(7,j)*P(7,j+1)+d(7,j)*P(6,j);dP(7,j)= P1(7,j)-P(7,j);endfor j=3:8P1(6,j)=c(6,j)*P(7,j)+a(6,j)*P(6,j-1)+e(6,j)*P(6,j)+b(6,j)*P(6,j+1)+d(6,j)*P(5,j);dP(6,j)= P1(6,j)-P(6,j);end%W2井底定压P1(6,9)=c(6,9)*P(7,9)+a(6,9)*P(6,8)+e(6,9)*P(6,9)+b(6,9)*P(6,10)+d(6,9)*P(5,9)-(3600e-9)*2*pi*K(6,9)*dt/u/C/Fai(6,9)/dx/dy/log(re/rw)*(P(6,9)-Pwf2);dP(6,9)= P1(6,9)-P(6,9);for j=3:7P1(5,j)=c(5,j)*P(6,j)+a(5,j)*P(5,j-1)+e(5,j)*P(5,j)+b(5,j)*P(5,j+1)+d(5,j)*P(4,j);dP(5,j)= P1(5,j)-P(5,j);endfor i=3:4for j=3:6%W1井定产if i==3&&j==5P1(i,j)=c(i,j)*P(i+1,j)+a(i,j)*P(i,j-1)+e(i,j)*P(i,j)+b(i,j)*P(i,j+1)+d(i,j)*P(i-1,j)-dt*Q1/24/C/Fai(i,j)/H(i,j)/dx/dy;dP(i,j)= P1(i,j)-P(i,j);elseP1(i,j)=c(i,j)*P(i+1,j)+a(i,j)*P(i,j-1)+e(i,j)*P(i,j)+b(i,j)*P(i,j+1)+d(i,j)*P(i-1,j);dP(i,j)= P1(i,j)-P(i,j);endendend%稳定条件压力变化最大值小于0.000001if max(max(abs(dP)))<0.000001break;endP=P1;tn=tn+1;end%将迭代时间步换算为时间天d=(tn-1)*dt/24;%计算W1井底流压。