数理方程基于matlab的数值解法
MATLAB方程与方程组的数值解ppt课件
4
第二节 MATLAB求解方程(组)的 函数及其用法
5
一、左除”\”与右除”/”
在MATLAB环境中，强烈建议使用左除”\”或者右除”/”解线性方程组 左除和右除是根据除号左侧还是右侧是分母而定的，方程系数矩阵 在未知数左侧，则用左除，反之用右除。使用左除”\”或者右除”/”的 好处是因为其对线性方程(组)的广泛适用性，当未知数个数大于方程 个数的时候，左除或右除会给出方程的特解，结合null函数，可以得 到通解。当未知数个数小于方程个数的时候，左除或右除会给出方 程的最小二乘解。
12
二、 MATLAB编程求解等额还款模型
➢ 给定月还款额、还款期数、贷款总额和利率计算到期剩余 贷款(AJfixPayment函数)
➢ 测试AJfixPayment 函数 ➢ 使用fsolve求出合适的月还款值，使得在120次还款后，
贷款余额为零。(SolveAJfixPayment函数)
13
r = roots(c)
其中输入参数： c: 多项式方程系数组成的行向量或者列向量，按降幂顺序排列。
函数输出参数： r: 多项式方程的解向量
参考： 【例11.2-4】8 Nhomakorabea四、 fsolve 函数
[x,fval,exitflag,output,jacobian] = fsolve(fun,x0,options) 其中输入参数：
Fun: 目标函数，一般用函数句柄形式给出 X0： 优化算法初始迭代解 Options: 参数设置（具体设置参考帮助文档） 输出参数： X: 最优解输出（或最后迭代解） Fval: 最优解（或最后迭代解）对应的函数值 Exitflag: 函数结束信息 （具体参考帮助文档 ） Output: 函数基本信息 包括迭代次数，目标函数最大计算次数，使 用的算法名称，计算规模 Jacobian：Jacobian矩阵（主要用来判断是否得到有效解）
基于MATLAB的常微分方程数值解法综述及经济模型
财经研究19基于MATLAB 的常微分方程数值解法综述及经济模型刘 欣 刘颖华 李海明摘 要：本文对常微分方程初值问题现有的数值解法进行了综述研究。
常用的数值解法：即欧拉法，改进欧拉法，龙格库塔方法，阿达姆斯公式等。
对较为典型的微分方程模型进行数值求解，探讨了数值算法在建模问题中的应用。
关键词：MATLAB 常微分方程 数值解法 模型 引言 在工程技术问题中，经常需要求解常微分方程的初值问题 ()()00,dy f x y dx y x y ⎧=⎪⎨⎪=⎩ （1） 微分方程的数值解法是一种离散化方法。
求这些值称为初值问题的一个数值解。
数值解法的基本思想是求初值问题(1)的解y （x）在一系列等距节点：0123n x x x x x <<<<<处的近似值：0123,,,,,n y y y y y 。
其中相邻两个节点间的距离1i i h x x +=-称为步长，即节点0 (0,1,2,)k x x kh k n =+= 。
1、数值求解法简介常微分方程的数值求解方法很多，常用的有欧拉(Euler)方法、休恩(Heun)方法、龙格－库塔(Runge-Kutta)法。
1.1 欧拉(Euler)方法根据泰勒(Taylor)定理，如果函数,,y y y '''连续，()()()0001,t y t hf t y y +=为欧拉逼近,将区间[a, b]分成M个等距的子区间，子区间的长度即为步长h=(b-a)/M，t k+1= t k +h。
重复上述欧拉逼近过程：()1112,y t hf y y +=；()2223,y t hf y y +=；这些点就是对曲线y=y(t)的逼近。
在(t 1, y 1)点，重复上述步骤，得到(t 2,y 2)，….，最后得到(t N , y N )。
步长取得越小，越逼近精确解。
1.2 休恩(Heun)方法对微分方程()y t f dt dy ,=从t 0到t 1积分，利用初始条件()00y t y =，得到 ()()()()()011010',t y t y dt t y dt t y t f tt t t -==⎰⎰→()()()()⎰+=1,01t t dt t y t f t y t y本方法需用欧拉方法作为预报值，并结合梯形公式(()()[]b f a f a b S +-=21)进行校正，即()()()1111,,2+++++=k k k k k k p t f y t f h y y 可见，用休恩方法要比用欧拉方法收敛快。
数值分析中求解线性方程组的MATLAB程序(6种)
数值分析中求解线性方程组的MATLAB程序（6种）1.回溯法（系数矩阵为上三角）function X=uptrbk(A,B)%求解方程组，首先化为上三角，再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));2.系数矩阵为下三角function x=matrix_down(A,b)%求解系数矩阵是下三角的方程组n=length(b);x=zeros(n,1);x(1)=b(1)/A(1,1);for k=2:1:nx(k)=(b(k)-A(k,1:k-1)*x(1:k-1))/A(k,k);end3.普通系数矩阵（先化为上三角，在用回溯法）function X=uptrbk(A,B)%求解方程组，首先化为上三角，再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));4.三角分解法function [X,L,U]=LU_matrix(A,B)%A是非奇异矩阵%AX=B化为LUX=B，L为下三角，U为上三角%程序中并没有真正解出L和U，全部存放在A中[N,N]=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);R=1:N;for p=1:N-1[max1,j]=max(abs(A(p:N,p)));C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=R(p);R(p)=R(j+p-1);R(j+p-1)=d;if A(p,p)==0'A is singular.No unique solution'break;endfor k=p+1:Nmult=A(k,p)/A(p,p);A(k,p)=mult;A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);endendY(1)=B(R(1));for k=2:NY(k)=B(R(k))-A(k,1:k-1)*Y(1:k-1);endX(N)=Y(N)/A(N,N);for k=N-1:-1:1X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);endL=tril(A,-1)+eye(N)U=triu(A)5.雅克比迭代法function X=jacobi(A,B,P,delta,max1);%雅克比迭代求解方程组N=length(B);for k=1:max1for j=1:NX(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)breakendendX=X';k6.盖斯迭代法function X=gseid(A,B,P,delta,max1);%盖斯算法，求解赋初值的微分方程N=length(B);for k=1:max1for j=1:Nif j==1X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1);elseif j==NX(N)=(B(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);elseX(j)=(B(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);endenderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)break;endendX=X';k。
用matlab求这个方程的数值解
用matlab求这个方程的数值解：y''-(1-y^2)y'+y=0, y(1)=2,y(2)=0, x?[0,20]，麻烦把程序编出来l84232525 11级分类：电脑/网络被浏览78次 2013.04.02用matlab求这个方程的数值解：y''-(1-y^2)y'+y=0, y(1)=2,y(2)=0,x?[0,20]，麻烦把程序编出来宝台木莲官采纳率：44% 10级 2013.04.03这是你们老师出的题目吧，这个例子在好多书上都出现过。
具体如下：function dydt = v dp1(t,y) dydt = [y(2); (1-y(1)^2)*y(2)-y(1)]; end 保存函数文件[t,y] = ode45(@vdp1,[0 20],[2; 0]); %This example uses @ to pass vdp1 as a function handle to ode45. T he resulting output is a column vector of time points t and a solution array y.Each row in y corresponds to a time returned in the corresponding row of t. 作图：You can simply use the plot command to view the solver output：plot(t,y(:,1),'-',t,y(:,2),' --') xlabel('time t'); ylabel('solution y'); legend('y_1','y_2')0.001y'''+0.11y''+y'+10y=10 初值X=0，y=0 要三阶微分方程的推导公式，如果有程序可以加分。
第八章基于MATLAB的科学计算—常微分方程数值解法
第八章基于MATLAB的科学计算—常微分方程数值解法第八章主要介绍了基于MATLAB的科学计算中的常微分方程数值解法。
常微分方程是数学中一个重要的领域，它研究的是一种未知函数的导数与
自变量之间的关系。
在实际问题中，常微分方程常常被用来描述物理系统、化学反应等现象。
在科学计算中，常微分方程的数值解法是一种常用的计算方式。
本章
首先介绍了常微分方程的基本概念和形式，包括一阶和高阶常微分方程的
定义和求解方法。
然后，详细介绍了数值解法中的欧拉方法、中点方法和
四阶龙格-库塔方法等。
欧拉方法是常微分方程数值解法中最简单的一种方法，它以初始条件
为基础，通过逐步逼近的方式求解方程。
中点方法在欧拉方法的基础上进
行了改进，通过使用两个点的平均斜率来逼近函数的斜率。
四阶龙格-库
塔方法是常微分方程数值解法中最常用的一种方法，它通过使用多个点的
斜率，以及加权平均值的方式来求解方程。
本章还介绍了MATLAB中的常微分方程求解工具箱，包括ode45、
ode23等函数，这些函数可以方便地求解各种不同类型的常微分方程。
MATLAB的求解工具箱提供了多种求解常微分方程的方法，用户可以根据
实际需求选择合适的方法。
总之，本章内容涵盖了基于MATLAB的科学计算中常微分方程数值解
法的基本概念和求解方法。
掌握这些知识，可以帮助科学计算中的研究人
员和工程师更好地理解和应用常微分方程数值解法，在实际问题中进行精
确的数值计算和模拟。
matlab化解方程
matlab化解方程用MATLAB求解方程是一种常见的数值计算方法，它可以帮助我们快速而准确地求解各种复杂的数学方程。
本文将介绍如何使用MATLAB 解决方程，并通过实例演示其应用。
我们要明确需要解决的方程是什么。
假设我们有一个一元二次方程，形如ax^2 + bx + c = 0，其中a、b、c是已知的系数。
我们的目标是求出方程的解x。
我们需要在MATLAB中定义方程的系数a、b、c。
我们可以使用变量来表示它们，例如a = 1，b = 2，c = 1。
接下来，我们可以使用MATLAB中的根函数来求解方程的解。
根函数可以接受一个多项式方程的系数作为输入，并返回方程的根。
在这个例子中，我们可以使用根函数来求解一元二次方程的解。
在MATLAB中，根函数的语法是roots(coef)，其中coef是一个包含方程系数的向量。
在我们的例子中，coef = [a b c]，即coef = [1 2 1]。
在MATLAB命令窗口中输入roots([1 2 1])，即可得到方程的解。
在这个例子中，方程的解为x = -1。
除了使用roots函数外，MATLAB还提供了其他一些求解方程的函数。
例如，如果我们需要求解非线性方程组，可以使用fsolve函数。
如果我们需要求解常微分方程，可以使用ode函数。
除了使用MATLAB内置的函数外，我们还可以使用MATLAB的符号计算工具箱来求解方程。
符号计算工具箱可以帮助我们进行符号计算，得到方程的精确解。
使用符号计算工具箱求解方程需要先定义方程的符号变量。
在MATLAB中，我们可以使用syms函数来定义符号变量。
例如，我们可以使用syms x来定义变量x为符号变量。
接下来，我们可以使用solve函数来求解方程。
solve函数可以接受一个或多个方程作为输入，并返回方程的解。
在我们的例子中，我们只有一个方程，即x^2 + 2x + 1 = 0。
我们可以使用solve函数来求解方程的解。
实验二_基于Matlab的微分方程数值解法
实验二微分方程数值解法一．实验原理及实验内容：对微分方程描述的控制系统，利用欧拉法、二阶龙格-库塔法、四阶龙格-库塔法分别编写M文件，进行数值计算和作图。
1．分别用欧拉法、二阶龙格-库塔法、四阶龙格-库塔法求下面系统的输出响应y(t)在0≤t≤1上，h=0.1时的数值解。
'2,(0)1=-=y y y要求保留4位小数，并将三种方法的结果与真解2=进行比较。
()ty t e-2.若为如2y y y==何编程计算？',(0)1二．实验仪器：计算机Matlab软件三．实验数据记录：程序一：disp('欧拉算法');y=1;h=0.1;for i=0:0.1:1disp(y);y=y+h*(-2*y);enddisp('欧拉算法');ydisp('精确解');yy=exp(-2*t)h=0.1;disp('函数的2阶数值解为');disp('y=');y=1;for t=0:h:1;disp(y);k1=-2*y;k2=-2*(y+k1*h);y(i+1)=y(i)+(k1+k2)*h*1/2;endh=0.1;disp('函数的4阶数值解为');disp('y=');y=1;for t=0:h:1;disp(y);k1=-2*y;k2=-2*(y+k1*h*1/2);k3=-2*(y+k2*h*1/2);k4=-2*(y+k3*h);y=y+h*1/6*(k1+2*k2+2*k3+k4); end>>程序2：t=0:0.1:1;n=length(t);y(1)=1;h=0.1;for i=1:n-1y(i+1)=y(i)+h*(y(i)*y(i)); enddisp('欧拉算法');ydisp('精确解');yy=exp(-2*t)h=0.1;disp('函数的2阶数值解为');disp('y=');y=1;for t=0:h:1;disp(y);k1=y*y;k2=(y+k1*h)^2;y=y+(k1+k2)*h*1/2;endh=0.1;disp('函数的4阶数值解为');disp('y=');y=1;disp(y);k1=y*y;k2=(y+k1*h*1/2)^2;k3=(y+k2*h*1/2)^2;k4=(y+k3*h)^2;y=y+h*1/6*(k1+2*k2+2*k3+k4); end。
matlab方程组数值解法
matlab方程组数值解法Matlab方程组数值解法随着科学技术的发展，数值计算在科学研究和工程实践中的应用越来越广泛。
对于复杂的数学模型，通过解析方法求得准确的解析解往往是困难的甚至不可能的。
因此，数值解法成为了求解这些问题的重要手段之一。
Matlab作为一种强大的数值计算工具，提供了多种数值解法来解决方程组的数值求解问题。
在Matlab中，求解方程组的数值解法主要包括直接法和迭代法两种。
直接法是指通过一系列直接计算来求解方程组的解，常见的方法有高斯消元法和LU分解法。
迭代法则是通过迭代计算来逼近方程组的解，常见的方法有雅可比迭代法和高斯-赛德尔迭代法。
高斯消元法是一种经典的直接法，它通过消元和回代的方式将方程组化为简化的三角方程组，然后通过回代计算得到解。
Matlab中提供了直接调用的函数，如"linsolve"函数，可以直接求解线性方程组。
对于非线性方程组，可以通过牛顿法等迭代法来求解。
LU分解法是另一种常用的直接法，它将方程组的系数矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积，然后通过回代计算得到解。
在Matlab中，可以使用"lu"函数进行LU分解，并通过"\"运算符求解线性方程组。
雅可比迭代法是一种简单而有效的迭代法，它通过迭代计算逐步逼近方程组的解。
在每一步迭代中，通过将方程组中的每个未知数的迭代解代入到方程组中的对应方程中，得到新的近似解。
通过多次迭代，可以得到逼近方程组解的解向量。
在Matlab中，可以使用"jacobi"函数进行雅可比迭代。
高斯-赛德尔迭代法是雅可比迭代法的改进版，它在每一步迭代中使用上一步迭代得到的未知数的新近似解。
这样可以更快地逼近方程组的解。
在Matlab中，可以使用"gauss_seidel"函数进行高斯-赛德尔迭代。
除了这些常见的数值解法外，Matlab还提供了其他一些数值求解函数，如"fsolve"函数可以求解非线性方程组，"ode45"函数可以求解常微分方程组等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数理方程数值解法与其在matlab软件上的实现张体强1026222 廖荣发1026226[摘要]数学物理方程的数值解在实际生活中越来越使用，首先基于偏微分数值解的思想上，通过matlab软件的功能，研究其数学物理方程的数值解，并通过对精确解和数值解进行对比，追究其数值解的可行性，在此，给出相关例子和程序代码，利于以后的再次研究和直接使用。
[关键字]偏微分方程数值解matlab 数学物理方程的可视化一：研究意义在我们解数学物理方程，理论上求数学物理方程的定解有着多种解法，但是有许多定解问题却不能严格求解，只能用数值方法求出满足实际需要的近似解。
而且实际问题往往很复杂，这时即便要解出精确解就很困难，有时甚至不可能，另一方面，在建立数学模型时，我们已作了很多近似，所以求出的精确解也知识推导出的数学问题的精确解，并非真正实际问题的精确解。
因此，我们有必要研究近似解法，只要使所求得的近似解与精确解之间的误差在规定的范围内，则仍能满足实际的需要，有限差分法和有限元法是两种最常用的求解数学物理方程的数值解法,而MATLAB 在这一方面具有超强的数学功能，可以用来求其解。
二：数值解法思想和步骤2.1:网格剖分为了用差分方法求解上述问题，将求解区域{}(,)|01,01x t x t Ω=≤≤≤≤作剖分。
将空间区间[0,1]作m 等分，将时间[0,1]区间作n 等分，并记1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。
分别称h 和τ为空间和时间步长。
用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。
2.2：差分格式的建立0u u t x∂∂-=∂∂………………………………（1） 设G 是,x t 平面任一有界域，据Green 公式(参考数学物理方程第五章)：()()Gu udxdt udt udx t xΓ∂∂-=--∂∂⎰⎰ 其中G Γ=∂。
于是可将（1）式写成积分守恒形式：()0udt udx Γ--=⎰ (2)我们先从（2）式出发构造熟知的Lax 格式设网格如下图所示图1取G 为以(1,)A j k +，(1,1)B j k ++，(1,1)C j k -+和(1,)D j k -为顶点的矩形。
T ABCD = A 为其边界，则()()()()()DABCABCDudt udx u dx u dx u dt u dt Γ--=-+-+-+-⎰⎰⎰⎰⎰ (3)右端第一个积分用梯形公式，第二个积分用中矩形公式，第三、四个积分用下矩形公式，则由（2）（3）式得到Lax-Friedrich 格式111111()202k k k k k j j j j j u u u u u hτ+-+-+-+-+=截断误差为()[]k k kj h j j R u L u Lu =-111111()22k k k k k k k j j j j j j j u u u u u u u h t xτ+-+-+-+-∂∂=+-+∂∂232223(),(0,0)26kkjj u u h O h j m k n t xττ∂∂=-=+≤≤≤≤∂∂ 所以Lax-Friedrich 格式的截断误差的阶式2()O h τ+ 令/s h τ=：则可得差分格式为111111(),(0,0)222k k k kk j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤ 0cos()(0)j j u x j m π=≤≤k+10cos(),cos(),(0)k k k m k u t u t k n ππ==-≤≤2.3差分格式的求解差分格式111111(),(0,0)222k k k k k j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤写成如下矩阵形式：1011121211110000022221100022221100000022220000000000000k k k k k k m m k m k m s s u u s s u u u us s u u +++---⎛⎫-+ ⎪⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪-+ ⎪ ⎪ ⎪⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-+ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭ ⎪ ⎪⎝⎭则需要通过对k 时间层进行矩阵作用求出k+1时间层。
对上面的矩阵形式我通过C++（或matlab ）编出如附录的程序求出数值解、真实解和误差。
例1：如下方程0,01,0 1.(,0)cos(),0 1.(0,)(1,)cos(),0 1.u ux t t xu x x x u t u t t t ππ∂∂-=≤≤≤≤∂∂=≤≤=-=≤≤ 利用 matlab 的数值解法结果并填入表中关系对比如下：表11/900,1/1000(0.9h s τ===）从上面可以看出，数值解精度相当高的三：matlab的在数学物理方程上简单的应用Matlab是一个强大的计算工具，超强的计算能力和绘图能力，下面几个例题来说明matlab数学物理上的应用例1：将函数1/（1－a）2在z=0 处展成幂级数。
解：>>syms a；>>Taylor (1/(1-a)^2,0)ans =1+2*a+3*a^2+4*a^3+5*a^4+6*a^5例2：写出函数f（x）＝1/(x2＋p2 )（a>0）的Fourier 变换式。
解：>>syms x w;>>syms a potitive>>f=1/(x^2+p^2);>>F=Fourier (f,x,w)F=pi*(signum(0,Re(p),0)*cosh(p*w)-2*Heaviside(w)*sinh(p*w)+sinh(p*w))/p例2：已知函数f（x）＝x3 e－x，试求取该函数的Laplace 变换,并对结果进行Laplace 反变换。
解：>>syms x w;>>f=x^3*exp(-x);>>F=laplace(f,x,w)F=6/(w+1)^4对得出的结果进行Laplace 反变换,从而有>>ilaplace(F)ans=x^3*exp(-x)利用手工方法对函数进行Fourier 变换和Laplace 变换，计算起来繁琐、复杂，且容易出错，利用MATLAB 快速、准确。
四：matlab解数学物理方程4.1：数值解法与精确解的可视化对比分析以下面的问题为例子（课本原题）根据上面的可建立方程如下：根据分离变量和差分其图形结果如下：图2 分离变量热传导图3 差分热传导其代码如下：图2x=0:0.1*pi:pi;y=0:0.4:10;[x,t]=meshgrid(x,y);u=0;m=length(j);%matlab可计算的最大数，相当于无穷for i=0:mu=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);图3u=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:20u(i,1)=i/20*pi;;end;for j=1:100u(1,j)=0;endfor j=1:99for i=2:19u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);endendfor j=1:100u(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(' 有限差分法解');从上面可以看出数值解法精度很高，图形基本完全一样的4.2 ：matlab实现数值解法以下面的方程为例基本步骤：(1)区域的离散或子区域的划分；(2)插值函数的选择；(3)方程组的建立；(4)方程组的求解。
a=input(' 请输入系数a 的值：');l=input(' 请输入长度l 的值：');M=input(' 请输入将区间[0,l]等分的个数M：');ot=input(' 请输入时间增量ot 的值：');n=input(' 请输入运行次数n 的值：');ox=l/M;x0=zeros(M+1,1);for ii=1:Mx0(ii+1)=ii*ox;endu=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-2B(ii)=1+2*r;A(ii)=-r;C(ii)=-r;S(ii)=u(ii+1,1);endB(M-1)=1+2*r;S(M-1)=u(M,1);u(1,2)=0;u(M+1, 2)=0;S(1,1)=S(1,1)+r*u(1,2);S(M-1,1)=S(M-1,1)+r*u (M+1,2);%追赶法S(1)=S(1)/B(1);T=B(1);k=2;while k~=MB(k-1)=C(k-1)/T;T=B(k)-A(k-1)*B(k-1);S(k)=(S(k)-A(k-1)*S(k-1))/T;k=k+1;endk=1;while k~=M-1S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k);k=k+1;endu(2:M,2)=S; %把结果放入矩阵u 中u(:,1)=u(:,2);% 过河拆桥end%计算精确值，存放在u 的第二列for x=0:Mu(x+1,2)=exp(-(pi*a/l)^2*n*ot)*sin(pi*x*ox/l); end%计算最大相对误差ez=zeros(M-1,1);for ii=2:Mez(ii-1)=abs(u(ii,1)-u(ii,2))/u(ii,2);endE=max(ez);fprintf (' 最后时刻数值解与精确解分别为:\n');disp (u);fprintf (' 差分法得到的结果与正确结果的最大相对误差为：');disp([num2str(E*100) '%']);%画二维图比较plot(x0,u(:,1),'r:',x0,u(:,2),'b-');legend(' 数值解',' 精确解')xlabel('x'),ylabel('u(x,t)')title(' 最后时刻热传导问题数值解与精确解比较')。