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

使用差分方法求解下面的热传导方程2(,)(,)(01,00.2,1)t xx T x t a T x t x t a =<<<<=初值条件：2(,0)44T x x x =-；边值条件：(0,)0(1,)0T t T t ==；使用差分公式1,,1,222(,)2(,)(,)2(,)()i j i j i j i j i j i jxx i j T x h t T x t T x h t T T T T x t O h h h -+--++-+=+≈,1,(,)(,)(,)()i j i j i j i jt i j T x t k T x t T T T x t O k kk++--=+≈上面两式带入原热传导方程,1,1,,1,22i j i ji j i j i jT T T T T kh +-+--+=令224kr h=，化简上式的,1,1,1,(12)()i j i j i j i j T r T r T T +-+=-++如下图：ix jt j编程MA TLAB 程序，运行结果如下1xtTfunction mypdesolution c=1;xspan=[0 1]; tspan=[0 0.2]; ngrid=[100 10]; f=@(x)4*x-4*x.^2; g1=@(t)0; g2=@(t)0;[T,x,t]=rechuandao(c,f,g1,g2,xspan,tspan,ngrid); [x,t]=meshgrid(x,t); mesh(x,t,T); xlabel('x') ylabel('t') zlabel('T')function [U,x,t]=rechuandao(c,f,g1,g2,xspan,tspan,ngrid) % 热传导方程：% Ut(x,t)=c^2*Uxx(x,t) a<x<b ts<t<tf % 初值条件： % u(x,0)=f(x)% 边值条件：% u(a,t)=g1(t)% u(b,t)=g2(t)%% 参数说明% c：方程中的系数% f：初值条件% g1,g2：边值条件% xspan=[a,b]：x的取值范围% tspan=[ts,tf]：t的取值范围% ngrid=[n,m]：网格数量，m为x网格点数量，n为t的网格点数量% U：方程的数值解% x,t：x和t的网格点n=ngrid(1);m=ngrid(2);h=range(xspan)/(m-1);x=linspace(xspan(1),xspan(2),m);k=range(tspan)/(n-1);t=linspace(tspan(1),tspan(2),n);r=c^2*k/h^2;if r>0.5error('为了保证算法的收敛，请增大步长h或减小步长k!')ends=1-2*r;U=zeros(ngrid);% 边界条件U(:,1)=g1(t);U(:,m)=g2(t);% 初值条件U(1,:)=f(x);% 差分计算for j=2:nfor i=2:m-1U(j,i)=s*U(j-1,i)+r*(U(j-1,i-1)+U(j-1,i+1));endend%%%%%%本文来自互联网%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%。

Fpg一维扩散方程の有限差分法——计算物理实验作业七 陈万物理学 2013 级 题目：编程求解一维扩散方程の解uD 2u 2 (0 x a 0 ,0 t t max ) tx u(x,t) |t 0 e x ua 1ub 1 nc 1(x 0)a 2ub 2 uc 2( x a 0 )n取 a 1 1,b 1 1, c 1 0, a 2 1, b 21, c 2 0, a 0 1.0,t max 。

输出 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 矩阵の对角元素dd = 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 矩阵の对角元素下边一行元素ee= -ones(1,n-1);e(1,n-1) = -b2;%求 A 矩阵の对角元素上边一行元素ff= -ones(1,n-1);f(1,1) = -b1;R = zeros(k,n);%求 R%追赶法求解for i = 2:kR(i,1) = (b1*P2-h*a1)*u(i -1,1)+b1*u(i -1,2)+2*h*c1;for j = 2:n-1R(i,j) = u(i -1,j-1)+2*P2*u(i -1,j)+u(i -1,j+1);endR(i,n) = b2*u(i -1,n-1)+( b2*P2-h*a2)*u(i -1,n)+2*h*c2;M = chase(e,d,f,R(i,:));u(i,:) = M';plot(x,u(i,:)); axis([0 a0 0 t_max]); pause(0.1)endoutput = u;%绘图比较解析解和有限差分解[X,T] = meshgrid(x,t);Z = exp(X+0.1*T);surf(X,T,Z),xlabel( 'x'),ylabel('t'),zlabel('u'),title( '解析解 '); figuresurf(X,T,u),xlabel('x'),ylabel('t'),zlabel('u'),title( '有限差分解 ');子程序 2：function M = chase(a,b,c,f)%追赶法求解三对角矩阵方程， Ax=f%a 是对角线下边一行の元素%b 是对角线元素%c 是对角线上边一行の元素%M 是求得の结果，以列向量形式保存n = length(b);beta = ones(1,n-1);y = ones(1,n);M = ones(n,1);for i = (n-1):(-1):1a(i+1) = a(i);end%将 a 矩阵和 n 对应beta(1) = c(1)/b(1);for i = 2:(n-1)beta(i) = c(i)/( b(i) -a(i)*beta(i -1) );endy(1) = f(1)/b(1);for i = 2:ny(i) = (f(i) -a(i)*y(i -1))/(b(i) - a(i)*beta(i-1));endM(n) = y(n);for i = (n-1):(-1):1M(i) = y(i) -beta(i)*M(i+1);endend结果：比较解析两图，结果令人满意。

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模拟昆明学院2015届毕业设计（论文）设计（论文）题目一维热传导问题的数值解法及其MATLAB模拟子课题题目无姓名伍有超学号 2所属系物理科学与技术系专业年级 2011级物理学2班指导教师王荣丽2015 年 5 月摘要本文介绍了利用分离变量法和有限差分法来求解一维传导问题的基本解，并对其物理意义进行了讨论。

从基本解可以看出，在温度平衡过程中，杠上各点均受初始状态的影响，而且基本解也满足归一化条件，表示在热传导过程中杆的总热量保持不变。

通过对一维杆热传导的分析，利用分离变量法和有限差分法对一维热传导进行求解，并用MATLAB 数学软件来对两种方法下的热传导过程进行模拟，通过对模拟所得三维图像进行取值分析，得出由分离变量法和有限差分法绘制的三维图基本相同，且均符合热传导过程中温度随时间、空间的变化规律，所以两种方法均可用来解决一维热传导过程中的温度变化问题。

关键词：一维热传导；分离变量法；有限差分法；数值计算；MATLAB 模拟AbstractIn this paper, the method of variable separation andfinite difference method are introduced to solve the problem of one-dimensional heat conduction problems, and the physical significance of numerical methods for heat conduction problems are discussed. From the basic solution, we can see the temperature on the bar are affected by the initial state during the process of temperature balance, and basic solution also satisfy the normalization condition which implied the invariance of the total heat in the bar during the heat conduction process. Through the analysis of the one-dimensional heat conduction, by taking use of variable separation method and finite difference method, we simulated the one-dimensional heat conduction problem by MATLAB. The three-dimensional images of the simulation results obtained by the method of separation of variables and finite difference method are similar to each other, and the temperature curve is in accordance with the law of temperature variation during heat conduction. Thus, we can go to the conclusion that both methods can be used to deal with the one-dimensional heat conduction problems.Keywords: One-dimensional heat conduction; method of variableseparation; finite difference method; numerical2method; MATLAB simulation目录第一章绪论11.1热传导的概念......................................................... .. (1)1.2热质的运动和传递......................................................... (1)第二章一维热传导问题的两种数值解法32.1一维热传导问题的初值问题32.2一维热传导问题的分离变量法42.3一维热传导问题的有限差分法63第三章一维有界杆热传导问题的MATLAB模拟9 3.1一维有界杆热传导问题93.2分离变量法的MATLAB模拟93.3有限差分法的MATLAB模拟12第四章总结与展望18参考文献19谢辞204第一章绪论1.1热传导的概念由于温度分布不均匀，热量从介质中温度高的地方流向温度低的地方称为热传导。

有限差分法的M a t l a b程序有限差分法的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-1x(i)=a+i*h;end % 定义网格点坐标for j=1:m-1y(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-2z=(-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);endu(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)>normnorm=abs(u(n-1,m-1)-z);endu(n-1,m-1)=z;% 对不靠近上下边界的网格点进行处理for j=m-2:-1:2% 对靠近左边界的网格点进行处理z=(-h^2*fun(x(1),y(j))+gun(a,y(j))+r*u(1,j+1)+r*u(1,j-1)+u(2,j))/s;if abs(u(1,j)-z)>normnorm=abs(u(1,j)-z);endu(1,j)=z;% 对不靠近左右边界的网格点进行处理for i=2:n-2z=(-h^2*fun(x(i),y(j))+u(i-1,j)+r*u(i,j+1)+r*u(i,j-1)+u(i+1,j))/s;if abs(u(i,j)-z)>normnorm=abs(u(i,j)-z);endu(i,j)=z;end% 对靠近右边界的网格点进行处理z=(-h^2*fun(x(n-1),y(j))+gun(b,y(j))+r*u(n-1,j+1)+r*u(n-1,j-1)+u(n-2,j))/s;if abs(u(n-1,j)-z)>normnorm=abs(u(n-1,j)-z);endu(n-1,j)=z;end% 对靠近下边界的网格点进行处理% 对左下角的网格点进行处理z=(-h^2*fun(x(1),y(1))+gun(a,y(1))+r*gun(x(1),c)+r*u(1,2)+u(2,1))/s;if abs(u(1,1)-z)>normnorm=abs(u(1,1)-z);endu(1,1)=z;% 对靠近下边界的除第一点和最后点外网格点进行处理for i=2:n-2z=(-h^2*fun(x(i),y(1))+r*gun(x(i),c)+r*u(i,2)+u(i+1,1)+u(i-1,1))/s;if abs(u(i,1)-z)>normnorm=abs(u(i,1)-z);endu(i,1)=z;end% 对右下角的网格点进行处理z=(-h^2*fun(x(n-1),y(1))+gun(b,y(1))+r*gun(x(n-1),c)+r*u(n-1,2)+u(n-2,1))/s;if abs(u(n-1,1)-z)>normnorm=abs(u(n-1,1)-z);endu(n-1,1)=z;% 结果输出if norm<=tolfid = fopen('FDresult.txt', 'wt');fprintf(fid,'\n********用有限差分法求解矩形域上Poisson方程的输出结果********\n\n'); fprintf(fid,'迭代次数: %d次\n\n',k);fprintf(fid,' x的值 y的值 u的值 u的真实值 |u-u(x,y)|\n');for i=1:n-1for j=1:m-1fprintf(fid, '%8.3f %8.3f %14.8f %14.8f %14.8f\n', [x(i),y(j),u(i,j),gun(x(i),y(j)),abs(u(i,j)-gun(x(i),y(j)))]);endendfclose(fid);break; % 用来结束while循环endk=k+1;endif k==N+1fid = fopen('FDresult.txt', 'wt');fprintf(fid,'超过最大迭代次数，求解失败！');fclose(fid);endclc[a1 a2 a3 a4] = textread('F:\aa.txt','%f %f %f %f');a = [a1 a2 a3];a=a';b=a4';[pa,mina,maxa,pb,minb,maxb]=premnmx(a,b);net =newrb(pa,pb,0,1.3,24,2);an =sim(net,pa);E = an - pb;m =sse(E)n = mse(E)[f1 f2 f3 f4]= textread('F:\bb.txt','%f %f %f %f');f = [f1 f2 f3];f=f';pf = tramnmx(f,mina,maxa);an2 = sim(net,pf);g =postmnmx(an2,minb,maxb);g= g';E2 = g- f4;mm =sse(E2)nn = mse(E2)。

一维热传导方程数值解法及matlab实现分离变量法和有限差分法一维热传导方程的Matlab解法：分离变量法和有限差分法。

问题描述：本实验旨在利用分离变量法和有限差分法解决热传导方程问题，并使用Matlab进行建模，构建图形，研究不同情况下采用何种方法从更深层次上理解热量分布与时间、空间分布关系。

实验原理：分离变量法：利用分离变量法，将热传导方程分解为两个方程，分别只包含变量x和变量t，然后将它们相乘并求和，得到一个无穷级数的解。

通过截取该级数的前n项，可以得到近似解。

有限差分法：利用有限差分法，将空间和时间分别离散化，将偏导数用差分代替，得到一个差分方程组。

通过迭代求解该方程组，可以得到近似解。

分离变量法实验：采用Matlab编写代码，利用分离变量法求解热传导方程。

首先设定x和t的范围，然后计算无穷级数的前n项，并将其绘制成三维图形。

代码如下：matlabx = 0:0.1*pi:pi;y = 0:0.04:1;x。

t] = meshgrid(x。

y);s = 0;m = length(j);for i = 1:ms = s + (200*(1-(-1)^i))/(i*pi)*(sin(i*x).*exp(-i^2*t));endsurf(x。

t。

s);xlabel('x')。

XXX('t')。

zlabel('T');title('分离变量法（无穷）');axis([0 pi 0 1 0 100]);得到的三维热传导图形如下：有限差分法实验：采用Matlab编写代码，利用有限差分法求解热传导方程。

首先初始化一个矩阵，用于存储时间t和变量x。

然后计算稳定性系数S，并根据边界条件和初始条件，迭代求解差分方程组，并将其绘制成三维图形。

代码如下：matlabu = zeros(10.25);s = (1/25)/(pi/10)^2;fprintf('稳定性系数S为:\n');disp(s);for i = 2:9u(i。

一维稳态导热问题数值模拟问题描述：设有一导热方程，022=+T dxTd ，边界条件为011dTx dx x T ⎧==⎪⎨⎪==⎩编写一段程序对此问题进行数值模拟。

解析：220d T T dx += 0011dT x dx x T ⎧==⎪⎨⎪==⎩ 1、用控制容积有限差分方法做出内部节点和边界节点的离散化方程：首先进行离散化，先确定节点，再确定控制容积。

将0-1划分为N 段，共N+1个节点，N 个控制容积，其中1x N∆=。

可以得到如下：对原方程建立差分方程，内部节点有：[()]0ew d dTT dx dx dx +=⎰ 0e wdT dT T x dxdx⇒-+∆=0P W E P P T T T T T x x x --⇒-+∆=∆∆1011P W E P P T T T T T N N N--⇒-+=1(2)1E W P E W P E W P P E E W W a a N N T NT NT N a a a a T a T a T N ==⎫⎧⇒-=+⎪⎪⇒⎬⎨=+-⎪⎪=+⎭⎩则转换为下式，：Pi i Ei Ei Wi Wi a T a T a T =+ i=2，….,N上式即为内部节点的离散化方程。

对于外部节点可有：1011i ii T T i T i N +==⎧⎨==+⎩综上可以得到内部节点和外部节点的离散化方程为：12111Pi i Ei Ei Wi Wi i i i a T a T a T i N T T i T i N +=+=⎧⎪==⎨⎪==+⎩，...,即为11(2)2111i Ei Wi i i i N T NT NT i N N T T i T i N +⎧-=+=⎪⎪⎨==⎪==+⎪⎩，...,上式不满足系数为负数，则可改用如下离散方程： 内部节点：*120E p P W P P T T T T T T N N N N----+= E w a a N == 1p E w a a a N =++ 12p a N N =+ *2p b T N=p p E E W W a T a T a T b =++ pi pi Ei Ei Wi Wi i a T a T a T b =++ *1112(2)()i i i P N T N T T T N N-++=++边界节点 1x= 11N T +=p p E E W W a T a T a T b =++E w a a N == 1p E w a a a N=++12p a N N=+*11112N N N N N N P a T a T a T T N ++--=++ *112(2)N N P N T N NT T N N-+=++边界节点 0x=0dTdx= (())0e P d dT T dx dx dx +=⎰ *1(2)012P E P P T T T T N N-+-=E a N = *1P b T N = 11++22P E a a N N N==p p E E a T a T b =+ *11221P a T a T T N =+ *1211()2P N T NT T N N-=+组成代数方程组：*12*11*111(+)1212(2)()212(2)1P i i i PN N P N T NT T i N N N T N T T T i N N N N T N NT T i N N N -+-⎧=+=⎪⎪⎪+=++≤≤⎨⎪⎪+=++=+⎪⎩写成矩阵方程组：*1*22*1*11+000021220001..0200 (2100202100)2P P N N P N P N N T N N T N N N T N T NN N N N T T T N N N NT N T N NN N N --⎡⎤-⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎡⎤-+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎢⎥+⎢⎥-+⎣⎦⎢⎥⎣⎦2、写出代数方程组的迭代求解程序: 用Matlab 编写如下求解程序；（1）高斯赛德尔迭代法（调用程序gauseidel 文件） function [x,n]=gauseidel(A,b,x0,eps,M) %高斯迭代格式 %线性方程组的系数：A%线性方程组中常数向量：b%迭代初始向量：x0%解的精度控制：eps%迭代步数控制：M%线性方程组的解：x%求出所需精度的解实际迭代步数：n if nargin==3eps=0.000001;M=10000;elseif nargin==4M=10000;endD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-L)\U;f=(D-L)\b;x=x0;n=0;tol=1;while tol>=epsx=G*x0+f;n=n+1;tol=norm(x-x0);x0=x;if (n>=M)disp ('Warning：’迭代次数过多，可能不收敛.') return;endend（2）主程序（demo文件）如下：N=input('请输入N值''\n')Tp=input('请输入Tp值''\n')x1=zeros(N,1)A0=zeros(N);A0(1,1)=N+1/(2*N);A0(1,2)=-N;A0(N,N-1)=-N;A0(N,N)=2*N+1/N;for i=2:N-1A0(i,i-1)=-N;A0(i,i)=2*N+1/N;A0(i,i+1)=-N;endb0=zeros(N,1);b0(1,1)=(1/N)*Tp;b0(N,1)=(2/N)*Tp+N;for i=2:N-1b0(i,1)=(2/N)*Tp;endA=A0;b=b0;x0=x1;[x,n]=gauseidel(A,b,x0)x=[x;1]t=(0:1/N:1)title('一维稳态导热问题空间温度分布图')xlabel('空间分布X')ylabel('温度分布T')hold onplot(t,x)T=1。

△t
n
nn
关于
t
的二阶中心差商[10]：
坠2u 坠x2
≈
uj+1
-2uj +uj-1 (△x)2
，对方
程进行离散。 离散后的方程为：
n n-1
n
nn
uj -uj △t
=a2
uj+1
-2uj +uj-1 (△x)2
。
令
：r=
a2·△t (△x)2
，即
n
n
n
n-1
（1＋2r)uj -r·uj+1 -r·uj-1 =uj 。 可化为矩阵形式：
摘 要：对于有界热传导齐次方程的混合问题，用分离变量法求解往往很复杂。 为了更好地
理解热传导方程的解，使用 MATLAB 软件将方程的解用图像表示出来。 通过区域转换的思想，
利用 MATLAB 编程实现一定区域内热传导方程的有限差分方法，数值表明了方法的可行性和
稳定性。
关键词：热传导方程；有限差分；MATLAB
方法， 把控制方程中的导数用网格节点上的函数值
的差商代替进行离散，从而 建立以网格节点上的值
为未知数的代数方程组。
1 求解热传导方程的基本思想
基本思想是把连续的定解区域用有限个离散点
构成的网格来代替， 这些离散点称作网格的节点；
把连续定解区域上的连续变量的函数用在网格上定
义的离散变量函数来近似； 把原方程和定解条件中
x0(ii+1)=ii*ox; end u=sin(pi*x0/l); % t=0 时 u(x,t)的值 r=a^2*ot/(ox)^2; for ii=1:n
%数据的输入 B=zeros(M-1,1);%存放系数矩阵主对角线元素 A=zeros (M-2,1);%存放系数矩阵主对角线元素下 方次对角线的元素 C=zeros (M-2,1);%存放系数矩阵主对角线元素上 方次对角线的元素 S=zeros(M-1,1);%存放右端的常数项 for ii=1:M-2

昆明学院2015届毕业设计（论文）设计（论文）题目一维热传导问题的数值解法及其MATLAB模拟子课题题目无姓名伍有超学号201117030225所属系物理科学与技术系专业年级2011级物理学2班指导教师王荣丽2015 年 5 月摘要本文介绍了利用分离变量法和有限差分法来求解一维传导问题的基本解，并对其物理意义进行了讨论。

从基本解可以看出，在温度平衡过程中，杠上各点均受初始状态的影响，而且基本解也满足归一化条件，表示在热传导过程中杆的总热量保持不变。

通过对一维杆热传导的分析，利用分离变量法和有限差分法对一维热传导进行求解，并用MATLAB 数学软件来对两种方法下的热传导过程进行模拟，通过对模拟所得三维图像进行取值分析，得出由分离变量法和有限差分法绘制的三维图基本相同，且均符合热传导过程中温度随时间、空间的变化规律，所以两种方法均可用来解决一维热传导过程中的温度变化问题。

关键词：一维热传导；分离变量法；有限差分法；数值计算；MATLAB 模拟AbstractIn this paper, the method of variable separation and finite difference method are introduced to solve the problem of one-dimensional heat conduction problems, and the physical significance of numerical methods for heat conduction problems are discussed. From the basic solution, we can see the temperature on the bar are affected by the initial state during the process of temperature balance, and basic solution also satisfy the normalization condition which implied the invariance of the total heat in the bar during the heat conduction process. Through the analysis of the one-dimensional heat conduction, by taking use of variable separation method and finite difference method, we simulated the one-dimensional heat conduction problem by MATLAB. The three-dimensional images of the simulation results obtained by the method of separation of variables and finite difference method are similar to each other, and the temperature curve is in accordance with the law of temperature variation during heat conduction. Thus, we can go to the conclusion that both methods can be used to deal with the one-dimensional heat conduction problems.Keywords: One-dimensional heat conduction; method of variable separation;finite difference method; numerical method; MATLAB simulation目录第一章绪论 (1)1.1热传导的概念 (1)1.2热质的运动和传递 (1)第二章一维热传导问题的两种数值解法 (3)2.1一维热传导问题的初值问题 (3)2.2一维热传导问题的分离变量法 (4)2.3一维热传导问题的有限差分法 (6)第三章一维有界杆热传导问题的MATLAB模拟 (9)3.1一维有界杆热传导问题 (9)3.2分离变量法的MATLAB模拟 (9)3.3有限差分法的MATLAB模拟 (12)第四章总结与展望 (18)参考文献 (19)谢辞 (20)第一章绪论1.1热传导的概念由于温度分布不均匀，热量从介质中温度高的地方流向温度低的地方称为热传导。

万方数据万方数据万方数据万方数据热传导方程有限差分法的MATLAB实现作者：史策作者单位：西安建筑科技大学,理学院,陕西,西安,710055刊名：咸阳师范学院学报英文刊名：JOURNAL OF XIANYANG NORMAL UNIVERSITY年，卷(期)：2009，24(4)被引用次数：0次1.曹钢,王桂珍,任晓荣.一维热传导方程的基本解[J].山东轻工业学院学报,2005,19(4):76-80.2.万正苏,方春华,张再云.关于热传导方程有限差分区域分解并行算法精度的注记[J].湖南理工学院学报(自然科学版),2007,20(3):12-14.3.StephenJ.Chapman.MATLAB编程[M].邢树军,郑碧波,译.北京:科学出版社,2008.4.田兵.用MATLAB解偏微分方程[J].阴山学刊,2006,20(4):12-13.5.王飞,裴永祥.有限差分方法的MATLAB编程[J].新疆师范大学学报(自然科学版),2003,22(4):21-27.6.王宝红.热传导方程的可视化探讨[J].忻州师范学院学报,2008,24(2):31-36.7.李先枝.热传导方程差分解法的最佳网格[J].河南大学学报(自然科学版),2004,34(3):16-18.8.赵德奎,刘勇.MATLAB在有限差分数值计算中的应用[J].四川理工学院学报,2005,18(4):61-64.9.谢焕田,吴艳.拉普拉斯有限差分法的MATLAB实现[J].四川理工学院学报,2008,21(3):1-2.10.南京大学数学系计算数学专业.偏微分方程数值解法[M].北京:科学出版社,1979.1.学位论文申卫东热传导方程有限差分区域分解算法研究2003区域分解算法是在并行机上求解偏微分方程数值解的一种较自然的方法.该方法先将偏微分方程求解区域划分为若干个子区域,然后在各个子区域并行求解.全文共五章.第一章为引言,简要介绍了热传导方程并行算法的概况及该文所讨论的基本内容.在第二章,我们在内边界点为等距分划的多子区域条件下,得到Dawson等人关于求解热传导方程区域分解算法差分解的误差估计.在第三章,我们以Saul'yev非对称格式作内边界处理,发展了新的区域分解算法,得到了差分解的先验误差估计,并与Dawson等人的算法作了比较.给出了关于算法计算精度的数值结果.在第四章,我们发展了一些新技术,在子区域的边界处采用小时间步长古典显式格式求解,构造了新的区域分解算法,得到了差分解的先验误差估计.给出了关于算法计算精度的数值结果.在第五章,我们在二维热传导方程求解上扩充了Dawson等人的区域分解算法.给出了关于算法计算精度的数值结果.第六章为该研究工作的主要结论.2.期刊论文张守慧.王文洽.ZHANG Shou-hui.WANG Wen-qia热传导方程有限差分逼近的数学Stencil及其新型迭代格式-山东大学学报（理学版）2006,41(6)将Stencil应用于偏微分方程有限元差分逼近过程,以两类差分格式为基础建立了求解热传导方程的两种新型迭代算法.此两种算法与经典的Jacobi方法同样具有并行的性质,但比Jacobi方法收敛快.给出的算例说明方法的适用性.3.期刊论文吕桂霞.马富明.Lü Guixia.Ma Fuming二维热传导方程有限差分区域分解算法-数值计算与计算机应用2006,27(2)本文讨论了一类数值求解二维热传导方程的并行差分格式.在这个算法中,通过引进内界点将求解区域分裂成若干子区域.在子区域间内界点上采用非对称格式计算,一旦这些点的值被计算出来,各子区域间的计算可完全并行.本文得到了稳定性条件和最大模误差估计.它表明我们的格式有令人满意的稳定性,并且有着较高的收敛阶.4.学位论文田源地下煤火三维数理模型正演数值模拟2006本文首先给出了几个地下煤火随空间、温度变化的动态和稳态热数学物理模型及其简化模型。

有限差分和pde 函数求解一维定态热传导方程分别用有限差分方法和pde 函数求解一维定态热传导方程，初始条件和边界条件，热扩散系数α=0.00001，22T T t x α∂∂=∂∂ （1） 求解过程：1. 用Tylaor 展开法推导出FTCS 格式的差分方程首先对T 进行泰勒展开得到如下两式子：23123123...2!3!23...2!3!nnnn n j j jj j nnnn n j jjjjttT T t x x TTx T T T t t t T T T x x x ++∆∆=+∆+++∆∆=+∆+++⎛⎫⎛⎫∂∂∂⎛⎫⎪ ⎪⎪∂∂∂⎝⎭⎝⎭⎝⎭⎛⎫⎛⎫∂∂∂⎛⎫ ⎪⎪ ⎪∂∂∂⎝⎭⎝⎭⎝⎭上述两个方程变换得：()1122323...23n nnn nn n j j j j jj j T T T T T t T t T o t t tt t t ++--⎛⎫⎛⎫∂∆∂∆∂⎛⎫=--=+∆ ⎪ ⎪ ⎪∂∆∂∂∆⎝⎭⎝⎭⎝⎭ （2）223123...23nnnn n j j jj j T T T x T x T x xx x --⎛⎫⎛⎫∂∆∂∆∂⎛⎫=-- ⎪ ⎪ ⎪∂∆∂∂⎝⎭⎝⎭⎝⎭()1232422342222...3!4!nnnn nnj j j jj j T T T T x T x T x x x x x x +-⎛⎫⎛⎫⎛⎫∂∂∆∂∆∂⎛⎫=--- ⎪ ⎪ ⎪ ⎪∂∆∆∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭ ()()21122222-n n n j j j T T T T o x x x+--+⎛⎫∂=+∆ ⎪∆∂⎝⎭ （3）将上述式子（2）（3）代入（1）得：12112222()nnn nn n n j j j j j j jT T T T T T T O t x t x t x αα+-+--+⎛⎫∂∂⎛⎫-=-+∆∆ ⎪ ⎪∂∂∆∆⎝⎭⎝⎭， （4）2. 方程的相容性和稳定性讨论：上述方程截项为：2233242334()...4...23!3!4!n n n n j jj j t T t T x T x TO t x t t x x α⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫∆∂∆∂∆∂∆∂ ⎪∆∆=--++ ⎪ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭， 由于(),0,0lim x t o t x ∆∆→∆∆=所以方程有相容性其经过傅里叶变换后，只有当,t x ∆∆满足下列条件时，方程具有较好的稳定性：2220sin 12m k xt x α∆∆≤≤∆ 其中m Nk Lπ=N 为节点个数，L 为边界长度由于：22sin sin 122m k x N x L π∆∆⎛⎫== ⎪⎝⎭所以当20.5ts x α∆=≤∆时方程具有稳定性3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。

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

下面是一个使用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');。

第五次作业（前三题写在作业纸上）一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.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 ，Txlabel('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<2fprintf('\n%s\n','有稳定性')elsefprintf('\n%s\n','没有稳定性')errorend%真实值计算[xe,te,Te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid);[xe,te]=meshgrid(xe,te);mesh(xe,te,Te);%画图，并且把坐标轴名称改为xe，te，Texlabel('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:jmaxrms=((T(j)-Te(j))^2/jmax)^(1/2)endfunction[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:nfor i=2:m-1for g=1:m-1Ue(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)) endendendfunction [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:nfor i=2:m -1U(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);endend5. 将温度随时间变化情况用曲线表示x t T6. 给出3000、9000、15000三个时刻的温度分布情况，对温度随时间变化规律做说明。

第四题完成
u(20,j)=u(19,j); end; disp(u); [x,t]=meshgrid(1:100,1:20); surf(x,t,u); xlabel('t'),ylabel('x'),zlabel('T'); title(' 有限差分法解'); 我们得到如图所示的热传导方程：
结论：
比较可得由以上两种方法作出的三维图形基本相同，符合热传导的热量分布 随时间和空间的变化规律
得到如图所示的热传导方程：
有限差分法
u=zeros(20,100); %t=1 x=pi 20 行 100 列 横坐标为 x 纵坐标为 t s=(1/100)/(pi/20)^2; fprintf(' 稳定性系数 S 为:\n'); disp(s); for i=1:20 u(i,1)=i/20*pi;; end; for j=1:100 u(1,j)=0; end for j=1:99 for i=2:19 u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); end end for j=1:100
模拟与仿真
根据课上所学知识，我们有如下方程：
u − a 2u = 0, 0 < x < l , t > 0 t xx u x x = l 0, t>0 = u x = 0 0, = 0< x<l = u t = 0 ϕ ( x ),
为便于解释做题，我们令： a=1 l=pi =x; 下面开始求解：
分离变量法 根据课上所讲
其中：
我们有如下代码：
x=0:0.1*pi:pi; y=0:0.4:10; [x,t]=meshgrid(x,y); u=0; m=length(j);%m u=u+8*(-1)^i/(pi*(2*i+1)^2)*(sin((2*i+1)/2*x).*exp(-(2*i+1)^2/4*t)); end; surf(x,t,u); xlabel('x'),ylabel('t'),zlabel('T'); title(' 分离变量法（无穷）'); disp(u);