1
科學計算(一)李義明老師2007年12月7日Heat Equation 程式使用說明Heat Equation 程式
(本程式使用向後差分法)
授課教師: 李義明教授
助教: 余紹銘、郭益廷
日期: 2007/12/7
2科學計算(一)李義明老師
2007年12月7日Heat Equation 程式使用說明function [u,x,t] = my_heat(t_0,t_f,x_0,x_f,M,N,type,T1,T2)% define the mesh in space
dx = (x_f-x_0)/N;
x = x_0:dx:x_f;
x = x';
% define the mesh in time
dt = (t_f-t_0)/M;
t = t_0:dt:t_f;
% define the ratio r
r = dt/dx^2;
fprintf('r = %g\n',r);
if r > 0.5
warning('r > 0.5')
dt = 0.499*dx^2;
t = t_0:dt:t_f;
fprintf('modified r = %g\n',dt/dx^2);
end
3科學計算(一)李義明老師
2007年12月7日Heat Equation 程式使用說明M = length(t);
N = length(x);
u = zeros(N,M);
if type == 0 % Dirichlet B.C
u(1,:) = T1;
u(N,:) = T2;
for i = 2:M % Time runs
A = zeros(N);
A(1,1) = 1;
u(1,i) = T1;
A(N,N) = 1;
u(N,i) = T2;
for j = 2:N-1
A(j,j) = 2*r + 1;
A(j,j+1) = -r;
A(j,j-1) = -r;
end
u(:,i) = A\u(:,i-1);
end
4科學計算(一)李義明老師
2007年12月7日Heat Equation 程式使用說明elseif type == 1 % Neumann at the left bounary, Dirichlet at the right boundary;u(N,:) = T2;
for i = 2:M % Time runs
A = zeros(N);
A(1,1) = 2*r+1;
A(1,2) = -2*r;
A(N,N) = 1;
u(N,i) = T2;
for j = 2:N-1
A(j,j) = 2*r + 1;
A(j,j+1) = -r;
A(j,j-1) = -r;
end
u(:,i) = A\u(:,i-1);
end
5科學計算(一)李義明老師
2007年12月7日Heat Equation 程式使用說明elseif type == 2 % Neumann at the right bounary, Dirichlet at the left boundary;u(1,:) = T1;
for i = 2:M % Time runs
A = zeros(N);
A(1,1) = 1;
u(1,i) = T1;
A(N,N) = 2*r+1;
A(N,N-1) = -2*r;
for j = 2:N-1
A(j,j) = 2*r + 1;
A(j,j+1) = -r;
A(j,j-1) = -r;
end
u(:,i) = A\u(:,i-1);
end
6科學計算(一)李義明老師
2007年12月7日Heat Equation 程式使用說明elseif type == 3 % Example 1: u(:,0) = sin(pi*x)
u(:,1) = sin(2*pi*x);
for i = 2:M % Time runs
A = zeros(N);
A(1,1) = 1;
u(1,i) = T1;
A(N,N) = 1;
u(N,i) = T2;
for j = 2:N-1
A(j,j) = 2*r + 1;
A(j,j+1) = -r;
A(j,j-1) = -r;
end
u(:,i) = A\u(:,i-1);
end
7科學計算(一)李義明老師
2007年12月7日Heat Equation 程式使用說明elseif type == 4 % Example 1: u(:,0) = sin(pi*x)
u(:,1) = zeros(N,1);
u(ceil(N/2),1) = 100;
for i = 2:M % Time runs
A = zeros(N);
A(1,1) = 1;
u(1,i) = T1;
A(N,N) = 1;
u(N,i) = T2;
for j = 2:N-1
A(j,j) = 2*r + 1;
A(j,j+1) = -r;
A(j,j-1) = -r;
end
u(:,i) = A\u(:,i-1);
end
8科學計算(一)李義明老師2007年12月7日Heat Equation 程式使用說明end
for i = 1:length(t)
figure(1), clf
plot(x,u(:,i),'.-');
axis([x(1) x(N) min(min(u)) max(max(u))])str = sprintf('t = %g sec',t(i));
xlabel('x'), ylabel('u(x)')
title(str)
pause(0.001);
end
9科學計算(一)李義明老師2007年12月7日Heat Equation 程式使用說明Numerical solution of the Black-
Scholes PDE
10
科學計算(一)李義明老師2007年12月7日Heat Equation 程式使用說明% EuPutExpl1.m
function price = EuPutExpl1(S0,X,r,T,sigma,Smax,dS,dt)% set up grid and adjust increments if necessary
M = round(Smax/dS);
dS = Smax/M;
N = round(T/dt);
dt = T/N;
matval = zeros(M+1,N+1);
vetS = linspace(0,Smax,M+1)';
veti = 0:N;
vetj = 0:M;
% set up boundary conditions
matval(:,N+1) = max(X-vetS,0);
matval(1,:) = X*exp(-r*dt*(N-veti));
matval(M+1,:) = 0;
% set up coefficients
a = 0.5*dt*(sigma^2*vetj -r).*vetj;
b = 1-dt*(sigma^2*vetj.^2 + r);
c = 0.5*dt*(sigma^2*vetj + r).*vetj;
11
科學計算(一)
李義明老師2007年12月7日Heat Equation 程式使用說明% solve backward in time
for i=N:-1:1
for j=2:M
matval(j,i) = a(j)*matval(j-1,i+1) + b(j)*matval(j,i+1)+ ...c(j)*matval(j+1,i+1);
end
end
% find closest point to S0 on the grid and return price % possibly with a linear interpolation
jdown = floor(S0/dS);
jup = ceil(S0/dS);
if jdown == jup
price = matval(jdown+1,1);
else
price = matval(jdown+1,1) + ...
(S0 -jdown*dS)*(matval(jup+1,1) -matval(jdown+1,1))/dS;end
12科學計算(一)李義明老師2007年12月7日Heat Equation 程式使用說明Example
>> [c,p] = blsprice(50,50,0.1,5/12,0.4);
>> p
p =
4.0760
>> EuPutExpl1(50,50,0.1,5/12,0.4,100,2,5/1200)ans =
4.0669
>> [c,p] = blsprice(50,50,0.1,5/12,0.3);
>> p
p =
2.8446
>> EuPutExpl1(50,50,0.1,5/12,0.3,100,2,5/1200)ans =
2.8288
有限差分法、有限单元法和有限体积法的简介 1.有限差分方法 有限差分方法(Finite Difference Method，FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。 2.有限元方法 有限元方法(Finite Element Method，FEM)的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。 有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。 在数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的
中心差分法的基本理论与程序设计 1程序设计的目的与意义 该程序通过用C语言（部分C++语言）编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。 2程序功能及特点 该程序采用C语言（部分C++语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。计算简便且在算法稳定的条件下，精度较高。 3中心差分法的基本理论 在动力学问题中，系统的有限元求解方程（运动方程）如下所示： ()()()() Ma t Ca t Ka t Q t ++= 式中，() a t分别是系统的结点加速度向 a t是系统结点位移向量，() a t和() 量和结点速度向量，,, M C K和() Q t分别是系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量，并分别由各自的单元矩阵和向量集成。 与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。 中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。通常的直接积分是基于两个概念，一是将在求解域0t T内的任何时刻t都应满足运动方程的要求，代之仅在一定条件下近似地满足运动方程，例如可以仅在相隔t?的离散的时间点满足运动方程；二是在一定数目的t?区域内，假设位移a、速度a、加速度a的函数形式。 中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为
求解线性方程组的直接解法 5.2LU分解 ① Gauss消去法实现了LU分解 顺序消元结束时的上三角矩阵U和所用的乘数,严格下三角矩阵。 将下三角矩阵的对角元改成1,记为L,则有A=LU, 这事实是一般的,我们不难从消去的第k个元素时的矩阵k行及k列元素的 历史得到这一点.因为从消元的历史有 u kj=a kj-m k1u1j- m k2u2j -…- m k,k-1u k-1,j, j=k,k+1,…,n m ik=(a ik-m i1u1k- m i2u2k -…-m i,k-1u k-1,k>/u kk i=k+1,k+2,…,n 于是a kj=m k1u1j+m k2u2j+…+m k,k-1u k-1,j+u kj, j=k,k+1,…,n a ik=m i1u1k+m i2u2k+…+m i,k-1u k-1,k+m ik u kk i=k+1,k+2,…,n 从前面两个式子我们可以直接计算L和U(见下段>.将矩阵分解为单位下 三角矩阵和上三角矩阵之积称为矩阵的LU分解.顺序消元实现了LU分 解,同时还求出了g, Lg=b的解. ②直接LU分解 上段我们得到(l ij=m ij> u kj=a kj-l k1u1j-l k2u2j -…- l k,k-1u k-1,j, j=k,k+1,…,n l ik=(a ik-l i1u1k-l i2u2k -…-l i,k-1u k-1,k>/u kk i=k+1,k+2,…,n 2 诸元素对应乘积,只不过算L的元素时还要除以同列对角元.这一规律很 容易记住.可写成算法(L和U可存放于A>: for k=1:n-1 for j=k:n u kj=a kj-l k1u1j-l k2u2j -…- l k,k-1u k-1,j end for i=k+1:n l ik=(a ik-l i1u1k-l i2u2k -…-l i,k-1u k-1,k>/u kk end end 这一算法也叫Gauss消去法的紧凑格式,可一次算得L,U的元素,不需逐步 计算存储.
有限元法与有限差分法的主要区别 有限差分方法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Ｔaylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观，表达简单，是发展较早且比较成熟的数值方法。对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑,可分为中心格式和逆风格式．考虑时间因子的影响,差分格式还可以分为显格式、隐格式、显隐交替格式等.目前常见的差分格式，主要是上述几种形式的组合,不同的组合构成不同的差分格式。差分方法主要适用于有结构网格,网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式:一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度,后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。有限元方法的基础是变分原理和加权余量法,其基本求解思想是把计算域划分为有限个互不重叠的单元,在每个单元内，选择一些合适的节点作为求解函数的插值点,将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式,借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解,整个计算域上总体的基函数可以看为由每个单元基函数组成的,则整个计算域内的解可以看作是由所有单元上的近似解构成．在河道数值模拟中,常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等.根据所采用的权函数和插值函数的不同,有限元方法也分为多种计算格式。从权函数的选择来说,有配置法、矩量法、最小二乘法和伽辽金法,从计算单元网格的形状来划分,有三角形网格、四边形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合同样构成不同的有限元计算格式。对于权函数，伽辽金(Gaｌeｒkiｎ）法是将权函数取为逼近函数中的基函数；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小;在配置法中，先在计算域内选取N个配置点。令近似解在选定的Ｎ个配置点上严格满足微分方程，即在配置点上令方程余量为０．插值函数一般由不同次幂的多项式组成,但也有采用三角函数或指数函数组成的乘积表示,但最常用的多项式插值函数。有限元插值函数分为两大类,一类只要求插值多项式本身在插值点取已知值,称为拉格朗日(Lａgranｇe)多项式插值;另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值,称为哈密特(Ｈｅrmite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状,一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元,常采用的插值函数为有Ｌa g ｒange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等．对于有限元方法，其基本思路和解题步骤可归纳为(１）建立积分方程,根据变分原理或方程余量与权函数正交化原理,建立与微分方程初边值问题等价的积分表达式,这是有限元法的出发点。(2）区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作,这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外,还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。(3）确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函
有限差分法 一、单变量函数： 用中心差分法（matlab程序见附录）计算结果如下： 图1 中心差分法
表1 数据对比 二、一维热传导： 在此取φ(x)=0，g1(t)= g2(t)=100-100*exp(-t)
问题描述： 已知厚度为l的无限大平板，初温0度，初始瞬间将其放于温度为100度的流体中，流体与板面间的表面传热系数为一常数。 试确定在非稳态过程中板内的温度分布。 (1)显式差分法： 图3 显式差分法 (2)隐式差分法： 图4 隐式差分法
小结：显式格式仅当时格式是稳定的。（其中称为网格比） 隐式格式从k层求k+1层时，需要求解一个阶方程组。而且隐式格式的稳定性对网格比没有要求，即为绝对稳定的。 三、Possion方程： 取f=1,R=1 图5差分法
图6 误差小结：观察误差曲面，其绝对误差数量级为
附Matlab程序： 第1题： %===========================Boundary Value Problem 1 clear;clc; A=[-2.01 1 0 0 0 0 0 0 0; 1 -2.01 1 0 0 0 0 0 0; 0 1 -2.01 1 0 0 0 0 0; 0 0 1 -2.01 1 0 0 0 0; 0 0 0 1 -2.01 1 0 0 0; 0 0 0 0 1 -2.01 1 0 0; 0 0 0 0 0 1 -2.01 1 0; 0 0 0 0 0 0 1 -2.01 1; 0 0 0 0 0 0 0 1 -2.01;]; c1=[0.1;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.9]; C=0.01*c1-1*[0;0;0;0;0;0;0;0;1]; y=A\C; x=0:0.1:1; yn=[0;y;1]; ye=2*(exp(x)-exp(-x))/(exp(1)-exp(-1))-x; figure(1); plot(x,yn,'*',x,ye); legend('numerical solution','exact solution') xlabel('x','fontsize',20); ylabel('y','fontsize',20); set(gca,'fontsize',18); figure(2); err=abs(ye'-yn); plot(x,err); legend('error') xlabel('x','fontsize',20); ylabel('y','fontsize',20); set(gca,'fontsize',18); 第2题： %========================Boundary Value Problem 1_Explicit %显式 clear;clc l=20;%板厚 h=1;%步长 J=l/h; T=50;%时间
利用有限差分法分析电磁场边界问题 在一个电磁系统中，电场和磁场的计算对于完成该系统的有效设计师极端重要的。例如，在系统中，用一种绝缘材料是导体相互隔离是，就要保证电场强度低于绝缘介质的击穿强度。在磁力开关中，所要求的磁场强弱，应能产生足够大的力来驱动开关。在发射系统中进行天线的有效设计时，关于天线周围介质中电磁场分布的知识显然有实质性的意义。 为了分析电磁场，我们可以从问题所涉及的数学公式入手。依据电磁系统的特性，拉普拉斯方程和泊松方程只能适合于描述静态和准静态（低频）运行条件下的情况。但是，在高频应用中，则必须在时域或频域中求解波动方程，以做到准确地预测电场和磁场，在任何情况下，满足边界条件的一个或多个偏微分方程的解，因此，计算电池系统内部和周围的电场和磁场都是必要的。 对电磁场理论而言，计算电磁场可以为其研究提供进行复杂的数值及解析运算的方法，手段和计算结果；而电磁场理论则为计算电磁场问题提供了电磁规律，数学方程，进而验证计算结果。常用的计算电磁场边值问题的方法主要有两大类，其每一类又包含若干种方法，第一类是解析法；第二类是数值法。对于那些具有最简单的边界条件和几何形状规则的（如矩形、圆形等）问题，可用分离变量法和镜像法求电磁场边值问题的解析解（精确解），但是在许多实际问题中往往由于边界条件过于复杂而无法求得解析解。在这种情况下，一般借助于数值法求解电磁场的数值解。 有限差分法，微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网络来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 差分运算的基本概念： 有限差分法是指用差分来近似取代微分，从而将微分方程离散成为差分方程组。于是求解边值问题即转换成为求解矩阵方程[5]。 对单元函数 ()x f而言，取变量x的一个增量x?=h，则函数()x f的增量可以表示为 ()x f? = ()h x f+-()x f 称为函数()x f 的差分或一阶差分。函数增量还经常表示为 ()x f? = ? ? ? ? ? + 2 h x f - ? ? ? ? ? - 2 h x f
追赶法求解三对角线性方程组 一 实验目的 利用编程方法实现追赶法求解三对角线性方程组。 二 实验内容 1、 学习和理解追赶法求解三对角线性方程组的原理及方法； 2、 利用MATLAB 编程实现追赶法； 3、 举例进行求解，并对结果进行分。 三 实验原理 设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵 11222=(1)(n 1)()()a c b a c A a n c b n a n ??????????--?????? ………… 这种方程组称为三对角线性方程组。显然，A 是上下半宽带都是1的带状矩阵。设A 的前n-1个顺序主子式都不为零，根据定理2.5的推论，A 有唯一的Crout 分解，并且是保留带宽的。 其中L 是下三角矩阵，U 是单位上三角矩阵。利用矩阵相乘法，可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ????????????????==?????--????????????……………
得到： 由上列各式可以得到L 和U 。 引入中间量y ，令 y Ux =，则有： 已知 L 和d ，可求得y 。 则可得到y 的求解表达式： 11/1 2,3,,()(1)*y()=()[()(1)]/y d l i n m i y i li i di y i di m i y i li ==-+=--… 1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui =*===≤≤+=+++≤≤-=?=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ?????????????????????????=??????---?????????????????? ……………
差分方程的解法分析及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 δ及其延迟表示,如果其值符合表达式则可合并处理.差分方程的完全解为
二阶常微分方程的中心差分求解 学校:中国石油大学（华东）理学院 姓名：张道德 一、 实验目的 1、 构造二阶常微分边值问题: 22,(),(), d u Lu qu f a x b dx u a u b αβ?=-+=<
11122 222222333222122112 100121012010012 00N N N u f q h h u f q h h h u f q h h h q u f h h ---???? ??+-???? ??? ???? ???????-+-? ?????? ???????????=-+? ?????? ???????????-???? ????????-+????? ?? ????? 可以看出系数矩阵为三对角矩阵，而对于系数矩阵为三对角矩阵的方程组可以用“追赶法”求解，则可以得出二阶常微分方程问题的数值解。 四、 举例求解 我们选取的二阶常微分方程边值问题为： 2 22242,01 (0)1,(1), x d u Lu x u e x dx u u e ?=-+=-<
《数值分析》课程设计追赶法解三对角方程组 院（系）名称信息工程学院 专业班级10普本信计 学号100111014 学生姓名刘银朋 指导教师张荣艳 2013 年05 月31日
数值分析课程设计评阅书 题目追赶法解三对角方程组 学生姓名刘银朋学号100111014 指导教师评语及成绩 指导教师签名： 年月日答辩评语及成绩 答辩教师签名： 年月日 教研室意见 总成绩： 教研室主任签名： 年月日
课程设计任务书 2012—2013学年第二学期 专业班级：10普本信息与计算科学学号：100111014 姓名：刘银朋 课程设计名称：数值分析Ⅰ、Ⅱ 设计题目：追赶法解三对角方程组 完成期限：自2013 年05月21 日至2013年05 月31日共10天 设计依据、要求及主要内容： 一、设计目的 理解追赶法,掌握追赶法的算法设计以及关于追赶法的分析和综合应用,能 够较熟练的应用Matlab软件编写求解追赶法的程序和应用Matlab软件数据库软 件. 二、设计内容 （1）认真挑选有代表性的三对角方程组. （2）认真梳理解三对角方程组的解题思路. （3）比较追赶法和高斯消去法的计算精度. 三、设计要求 1．先用Matlab数据库中的相应的函数对选定的方程,求出具有一定精度的解. 2．然后使用所用的方法编写Matlab程序求解. 3.对于使用多个方程解同意问题的,在界面上要设计成菜单的形式. 计划答辩时间：2013年06 月 5 日 工作任务鱼工作量要求： 查阅文献资料不少于3篇,课程设计报告1篇不少于3000字. 指导教师（签字）：教研室主任（签字）： 批准日期：2013 年05 月20 日
第三章 有限差分法 函数()f x ，x 为定义在区间[]a b ,上的连续变 量。将区间[]a b ,等分成n 份，令()h b a n =-称为 步长，x 在这些离散点处的取值为 x a ih i =+ ()i n =01,,, 称为节点。函数()f x 在这些节点处的差值 ()()()()()() f x h f x f x f x h f x h f x h i i i i i i +---+--??? ?? (5-1) 分别称为一阶向前、向后和中心差分，可以用它 们作为函数()f x 在x i 处的微分近似值。这些差分 与相应x 区间的比值 ()()[] ()() [ ] ()()[] 1 1 1 2h f x h f x h f x f x h h f x h f x h i i i i i i +---+--?????? ??? (5-2) 分别称为一阶向前、向后和中心差商，可以用它 们作为函数()f x 在x i 处的导数近似值。完全类似 地可以定义高阶差商，例如常用的二阶中心差商 ()()()[] 1 22h f x h f x f x h i i i +-+- (5-3) 可以作为函数()f x 在x i 处的二阶导数近似值。 §3.1 常微分方程初值问题的差分解法 考虑电学中的一个问题：如图5-1。研究 电容器上的电荷随时间的变化规律。 图5-1 RC 放电回路 这个问题对应的微分方程及其定解条件为：
d d Q t Q RC Q Q t =-=??? ??=00 (5-4) 这是一阶微分方程的初值问题，它的解析解为 Q Q e t RC =-0 (5-5) 一、欧拉(Euler )折线法 求解下列普遍形式的一阶微分方程的初值 问题： ()[]()'=∈=?????y f x y x a b y a y ,,0 (5-6) 首先，将区间[]a b ,等分n 份，取值 a x x x b n =<<<=01 ，步长h x x i i =-+1。 然后，用一阶向前差商近似一阶导数，即 ()() ()()[] y x y x h y x f x y x i i i i i +-≈'=1, (5-7) 简记()y x y i i ≈，则式(5-7)可以写成差分格式： ()y y h f x y i i i i +=+?1, ()i n =-011,,, (5-8) 此即向前欧拉差分格式。这是一个递推计算格式， 从区间左端点即式(5-6)中的初始条件出发，按式 (5-8)依次可以算到区间右端点，得到的 y y y n 12,,, 就是原方程解()y x 的近似值。 应用式(5-8)计算RC 放电方程(5-4)，按SI 单 位制，取Q 010=，RC =8，时间步长h =1，计 算结果如下：
中心差分法计算单自由度体系的动力反应 北京工业大学结构工程 组员:胡建华 S201204111 马 恒 S201204112 陈相家S201204083 张力嘉S201204022 0前言 时域逐步积分法是数值分析方法，它只假设结构本构关系在一个微小的时间步距内是线性的。时域逐步积分法是结构动力问题中一个得到广泛研究的课题，并在结构动力反应计算中得到广泛应用。由于引进的假设条件不同，可以有各种不同的方法，比如中心差分法，线性加速度法，平均常加速度法，Wilson-θ法等，其中中心差分法精度最高。在本文中，通过编制中心差分法计算单自由度体系的动力反应的程序来了解其应用及稳定性。 1中心差分法原理 中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。 中心差分法只在相隔t ?一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ?=?，则速度与加速度的中心差分近似为： t u u u i i ?-= -+? 21 1 (a) 2 1 12t u u u u i i i ?+-= -+? ? (b) 而离散时间点的运动为 )(),(),(i i i i i i t u u t u u t u u ? ?? ?? ? === ( =i 0,1,2,3,……) 由体系运动方程为：0)()()(=++? ? ?t ku t u c t u m i (c) 将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t 时刻的运动方程： 0221 12 11=+?-+?+--+-+i i i i i i ku t u u c t u u u m （d ）
3）三对角形线性方程组 123456789104100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014x x x x x x x x x x -????????--????????--????--????????--????--????????--????--???????--??????-???? 7513261214455????????-?? ?? ??=??-?? ???? -?? ?????? ???-?? *(2,1,3,0,1,2,3,0,1,1)T x =--- 二、数学原理 设系数矩阵为三对角矩阵 1 122233111000000000 000000 n n n n n b c a b c a b A a b c a b ---?? ? ? ?= ? ? ? ? ?? ? 则方程组Ax=f 称为三对角方程组。 设矩阵A 非奇异，A 有Crout 分解A=LU ，其中L 为下三角矩阵，U 为单位上三角矩阵，记 1 122 233 1 10 00010 000 0001000 000100,00000000 00 0001n n n n b L U γαβγββγβ--???? ? ? ? ? ? ??== ? ? ? ? ? ? ? ? ? ??? ? ? ? 可先依次求出L ，U 中的元素后，令Ux=y ，先求解下三角方程组Ly=f 得出y ，再求解上三角方程组Ux=y 。
事实上，求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起，使算法更为紧凑。其计算公式为： 1111, 1111 ,111 ,2,3,,,1,2,,1i i i i i i i i i i i i i i n n i i i i c f b y i n c a b a f y y x y i n n x y x βγββαβγγβαβγ--+? ===?? =?? ?==-= ??? -?=?? =??=--?=-??对对(*) 三、程序设计 function x=chase(a,b,c,f) %求解线性方程组Ax=f,其中A 是三对角阵 %a 是矩阵A 的下对角线元素a(1)=0 %b 是矩阵A 的对角线元素 %c 是矩阵A 的上对角线元素c(n)=0 %f 是方程组的右端向量 n=length(f); x=zeros(1,n);y=zeros(1,n); d=zeros(1,n);u= zeros(1,n); %预处理 d(1)=b(1); for i=1:n-1 u(i)=c(i)/d(i); d(i+1)=b(i+1)-a(i+1)*u(i); end %追的过程 y(1)=f(1)/d(1); for i=2:n y(i)=(f(i)-a(i)*y(i-1))/d(i); end %赶的过程 x(n)=y(n); for i=n-1:-1:1 x(i)=y(i)-u(i)*x(i+1); end
例题：已知差分方程51 (2)(1)()(+1)+0.5()66 x k x k x k r k r k +-++=，其中r (k )=1，k ≥0，x (0)=1， x (1)=2。 (1) 试由迭代法求其全解的前5项； (2) 分别由古典法求其零输入解、零状态解，以及全解； (3) 用Z 变换法求解差分方程。 解：注：解题过程中出现的下标“zi ”和“zs ”分别表示零输入条件和零状态条件。 1. 迭代法 题目中给出的条件仅仅是零输入初始条件，进行迭代求解时的初始条件应该是全解初始条件。 (1) 零输入初始条件 本题已给出零输入时的两个初始条件x zi (0)=1，x zi (1)=2。 (2) 零状态初始条件 取k =-2时，则51 (0)(1)(2)(1)0.5(2)66x x x r r --+-=-+-，得x zs (0)=0； 取k =-1 时，则51 (1)(0)(1)(0)0.5(1)66 x x x r r -+-=+-，求得x zs (1)=1。 (3) 全解初始条件 x (0)= x zi (0)+ x zs (0)=1； x (1)= x zi (1)+ x zs (1)=3。 (4) 根据求出的全解x (0)和x (1)，利用迭代法求解 取k =0时，则51(2)(1)(0)(1)0.5(0)66x x x r r -+=+，求得23(2)6x =； 取k =1时，则51(3)(2)(1)(2)0.5(1)66x x x r r -+=+，求得151 (3)36x =； 取k =2时，则51(4)(3)(2)(3)0.5(2)66x x x r r -+=+，求得941 (4)216 x =。 2. 古典法 (1) 零输入解 令输入为零，则得齐次方程 51 (2)(1)()066 x k x k x k +-++= (a) 根据差分方程定义的算子()()n d x k x k n =+，可得它的特征方程251 066 d d -+= 求得特征根为： 112d = ，21 3 d =
中心差分法计算程序编程 姓名：张泽伟 学号： 电话： 一、中心差分法程序原理说明 1.1 中心差分法思路 中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。 1.2 中心差分法原理 中心差分法只在相隔t ?一些离散的时间区间内满足运动方程，其基于有限 差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ?=?， 则速度与加速度的中心差分近似为： t u u u i i ?-= -+?211 (a) 21 12t u u u u i i i ?+-= -+?? (b) 而离散时间点的运动为 ) (),(),(i i i i i i t u u t u u t u u ??????=== ( =i 0,1,2,3,……) 由体系运动方程为：0 )()()(=++???t ku t u c t u m i (c) 将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t 时 刻的运动方程： 02211211=+?-+?+--+-+i i i i i i ku t u u c t u u u m （d ） 在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t 以前时刻的运动已知，则 可以把已知项移到方程的右边，整理得到： 12212)2()2()2(-+?-?-?--=?+?i i i u t c t m u t m k u t c t m (e)
由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t 时刻的运动，如果需 要可以用式（a ）和式（b ）求得体系的速度和加速度。 1.3 初始条件转化 假设给定的初始条件为 ),0(), 0(00? ?==u u u u （g ） 由式（g ）确定1-u 。在零时刻速度和加速度的中心差分公式为： t u u u ?-= -?2110 （h ） ` 21 0102t u u u u ?+-=-?? （i ） 将式（i ）消去1u 得：020012???-?+?-=u t u t u u （j ） 而零时刻的加速度值0??u 可以用t ＝0时的运动方程 0000=++???ku u c u m 确定 即 )(1000ku u c m u --=?? ? （k ） 这样就可以根据初始条件 00,?u u 和初始荷载0P ，就可以根据上式确定1-u 的 值。 1.4 中心差分法编程思路 ① 基本数据准备和初始条件计算： )(1000ku u c m u --=?? ? 020012??? -?+?-=u t u t u u ② 计算等效刚度和中心差分计算公式中的相关系数： t c t m k ?+?=22
题目：Newton-Raphson 迭代法 （1）计算原理 （2）编出计算机程序 （3）给出算例（任意题型） (1)计算原理： 牛顿-拉夫森(Newton-Raphson)迭代法也称为牛顿迭代法，它是数值分析中最重要的方法之一，它不仅适用于方程或方程组的求解，还常用于微分方程和积分方程求解。 用迭代法解非线性方程时，如何构造迭代函数是非常重要的，那么怎样构造的迭代函数才能保证迭代法收敛呢？牛顿迭代法就是常用的方法之一，其迭代格式的来源大概有以下几种方式： 1设()[]2,f x C a b ∈，对()f x 在点[]0,x a b ∈，作泰勒展开： 略去二次项，得到()f x 的线性近似式：()()()()000f x f x f x x x '≈+- 由此得到方程()0f x =的近似根(假定()00f x '≠)，() () 000f x x x f x =-' 即可构造出迭代格式(假定()00f x '≠)：() () 1k k k k f x x x f x +=- ' 这就是牛顿迭代公式，若得到的序列{}k x 收敛于α，则α就是非线性方程的根。 2 牛顿迭代法 牛顿切线法，这是由于()f x 的线性化近似函数()()()()000l x f x f x x x '≈+-是曲线()y f x =过点()()00,x f x 的切线而得名的，求()f x 的零点代之以求() l x !2))((''))((')()(2 0000x x f x x x f x f x f -+ -+= ξ
的零点，即切线与x 轴交点的横坐标，如左图所示，这就是牛顿切线法的几何解释。实际上，牛顿迭代法也可以从几何意义上推出。利用牛顿迭代公式，由 k x 得到1k x +，从几何图形上看，就是过点()(),k k x f x 作函数()f x 的切线k l ，切线k l 与x 轴的交点就是1k x +，所以有()() 1 k k k k f x f x x x +'=-，整理后也能得出牛顿迭 代公式： 3 要保证迭代法收敛，不管非线性方程()0f x =的形式如何，总可以构造： 作为方程求解的迭代函数。因为： 而且 在根附近越小，其局部收敛速度越快，故可令： 若0(即根不是0的重根)，则由得： ， 因此可令 ，则也可以得出迭代公式： 。 4 迭代法的基本思想是将方程改写成等价的迭代形式，但随之而来的问题却是迭代公式不一定收敛，或者收敛的速度较慢。运用前述加速技巧，对于简单迭代过程 ，其加速公式具有形式： ，其中 记，上面两式可以合并写成： 这种迭代公式称作简单的牛顿公式，其相应的迭代函数是： 。 需要注意的是，由于是的估计值，若取，则实际上便是的估计值。假设，则可以用代替上式中的， 就可得到牛顿法的迭代公式： 。 )(')(1k k k k x f x f x x - =+)()()(x f x k x x x -==?)0)((≠x k )(')()()('1)('x f x k x f x k x --=?) ('x ?α0)('=α?≠)('αf α=)(x f 0)('=α?)('1 )(ααf k = )('1 )(x f x k = )(')(1k k k k x f x f x x - =+0)(=x f )(x x ?=)(1n n n x f x x +=+θθ?--= +1)(1n n n x x x ) (111n n n x x x --+=++θθ )(1 n n x x ?=+1-=θL L x f x x n n n )(1- =+L x f x x )()(- =?L )('x ?)()(x f x x +=?)('x ?)('x f 0)('≠x f )('x f L )(')(1n n n n x f x f x x - =+
中心差分法求解单自由度体系的自由振动问题 前言 时域逐步积分法是根据运动方程，引进某些假设，建立由t 时刻状态向量i u 、i u ?、i u ? ?到t ＋t ?时刻的状态向量1+i u 、1+?i u 、1+??i u 的递推关系，从而从t ＝0时刻的初始状态向量0u 、0?u 、0? ?u 出发，逐步求出各时刻的状态向量，由于引进的假设条件不同，可以有各种不同的方法，下面主要介绍一种时域逐步积分方法－中心差分法。 中心差分法(central difference method)原理[1] 中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。 中心差分法只在相隔t ?一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ?=?，则速度与加速度的中心差分近似为： t u u u i i ?-=-+?211 (a) 2 112t u u u u i i i ?+-=-+?? (b) 而离散时间点的运动为 )(),(),(i i i i i i t u u t u u t u u ??????=== ( =i 0,1,2,3,……) 由体系运动方程为：0)()()(=++???t ku t u c t u m i (c) 将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t 时刻的运动方程： 02211211=+?-+?+--+-+i i i i i i ku t u u c t u u u m （d ） 在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到： 12212)2()2()2(-+?-?-?--=?+?i i i u t c t m u t m k u t c t m (e)