下载文档原格式
偏微分方程的数值解法上机习题七题目：求解边值问题：()()()()()()22sin cos cos sin ,0,10,10,,.x y u e x y x y x y G u x y G ππππππ+⎧∆=+⎪⎪∈=⨯⎨⎪=∈∂⎪⎩，，取步长1/64,1/128,h k == 作五点差分格式。
用Jacobi 迭代，Gauss-Seidel 迭代和SOR 迭代，（取opt ϖϖ= ）求解差分方程，以前后两次重合到小数点后四位的迭代值作为解的近似。
比较三种解法的迭代次数以及差分解()(),1/64,1/128h u x y h = 与精确解()(),sin sin x y u x y e x y πππ+= 的精度。
求解：1 给定步长，建立方程组系数矩阵和常数项。
2 计算出迭代法的迭代矩阵、常数项及最佳松弛因子3 进行迭代计算，输出迭代次数4 求解精确解并作出精确解和迭代解对应图像5 作出误差图像程序：clc;clear;n=input('请输入等分数 n=');h=1/n;%---------------------建立常数项b--------------------------for i=1:(n-1)for j=1:(n-1)b(j+(i-1)*(n-1),1)=2*pi^2*exp(pi*(i*h+j*h))*(sin(pi*i*h)*cos(pi*j*h)+cos(pi*i* h)*sin(pi*j*h));endend%---------------------建立系数矩阵A--------------------------j=0;for i=1:(n-1)^2if i>=nj=j+1;A(i,j)=1/h^2;A(j,i)=1/h^2;endA(i,i)=-4/h^2;endj=0;for i=2:(n-1)^2j=j+1;if mod(j,n-1)~=0A(j,i)=1/h^2;endend%--------------迭代法矩阵----------------D=diag(diag(A));L=tril(A,-1);R=triu(A,1);N=20000;%迭代最大次数%------------Jacboi迭代法迭代矩阵及常数项-----------------G=-inv(D)*(L+R); H=inv(D)*b;%------------Guass-Seidel迭代法迭代矩阵及常数项------------M=inv(D+L);R1=-M*R;M=M*b;%-----------SOR超松弛迭代法迭代矩阵、常数项及最佳松弛因子-------------I=speye((n-1)^2);B=I-inv(D)*A;q=eig(B);q=q(1);w=2/(1+sqrt(1-q^2));P=inv(D-w*(-L));Aw=P*(w*(-R)+(1-w)*D);P=w*P*b;%------------------------迭代法计算并作图-----------------------------str={'1.Jacobi迭代的迭代次数是','2.Guass-Seidel迭代的迭代次数是','3.SOR超松弛法的迭代次数是'};for k=1:3U=speye((n-1)^2,1);U0=zeros((n-1)^2,1);U0=sparse(U0);time=0;while norm(U-U0)>1e-4U0=U;switch kcase 1U=G*U0+H;case 2U=R1*U0+M;case 3U=Aw*U0+P;endtime=time+1;if i>Nbreakendenddisp([char(str{k}),num2str(time)]);for a=1:n-1for b=1:n-1switch kcase 1Z_1(a,b)=U((n-1)*(a-1)+b,1);case 2Z_2(a,b)=U((n-1)*(a-1)+b,1);case 3Z_3(a,b)=U((n-1)*(a-1)+b,1);endendendend%-----------------------作方程解图需要---------------------------------x=1/n:1/n:(n-1)/n;y=x;[X,Y]=meshgrid(x,y);figure(1);%-------------------------精确解计算并作图----------------------------- Z=exp(pi*(X+Y)).*sin(pi*X).*sin(pi*Y);subplot(2,2,1);mesh(X,Y,Z);title('精确解图像');xlabel('x');ylabel('y');zlabel('u');%-----------------------迭代法图像-------------------------------- subplot(2,2,2);mesh(X,Y,Z_1);title('Jacobi迭代图像');xlabel('x');ylabel('y');zlabel('u');subplot(2,2,3);mesh(X,Y,Z_2);title('Guass-Seidel迭代图像');xlabel('x');ylabel('y');zlabel('u'); subplot(2,2,4);mesh(X,Y,Z_3);title('SOR超松弛法图像');xlabel('x'); ylabel('y'); zlabel('u');%------------------------迭代法误差图像---------------------------- figure(2);subplot(2,2,1);mesh(X,Y,Z_1-Z);title('Jacobi迭代误差图像');xlabel('x');ylabel('y');zlabel('u'); subplot(2,2,2);mesh(X,Y,Z_2-Z);title('Guass-Seidel迭代误差图像');xlabel('x');ylabel('y');zlabel('u'); subplot(2,2,3);mesh(X,Y,Z_3-Z);title('SOR超松弛误差图像');xlabel('x'); ylabel('y'); zlabel('u');运行结果（1） 步长为1/64a) 迭代次数请输入等分数 n=641.Jacobi 迭代的迭代次数是382.Guass-Seidel 迭代的迭代次数是283.SOR 超松弛法的迭代次数是26b) 解的图像x精确解图像yuxJacobi 迭代图像yuxGuass-Seidel 迭代图像yuxSOR 超松弛法图像yuc) 误差图像xJacobi 迭代误差图像yuxGuass-Seidel 迭代误差图像yuxSOR 超松弛误差图像yu（2） 步长为1/32a) 迭代次数请输入等分数 n=321.Jacobi 迭代的迭代次数是392.Guass-Seidel 迭代的迭代次数是293.SOR 超松弛法的迭代次数是27b) 解的图像x精确解图像yuxJacobi 迭代图像yuxGuass-Seidel 迭代图像yuxSOR 超松弛法图像yuc) 误差图像xJacobi 迭代误差图像yuxGuass-Seidel 迭代误差图像yuxSOR 超松弛误差图像yu（3） 步长为1/16a) 迭代次数请输入等分数 n=161.Jacobi 迭代的迭代次数是382.Guass-Seidel 迭代的迭代次数是283.SOR 超松弛法的迭代次数是26b) 解的图像x精确解图像yuxJacobi 迭代图像yuxGuass-Seidel 迭代图像yuxSOR 超松弛法图像yuc) 误差图像xJacobi 迭代误差图像yuxGuass-Seidel 迭代误差图像yuxSOR 超松弛误差图像yu结果分析1从迭代次数来看，Jacobi 迭代法最高，Guass-Seidel 迭代法次之，SOR 超松弛法最少。
L试讨论逼近对蘇程詈+若。
的差分沁1)2)q1 二：行口匚1）解：设点为（X ? ，/曲）屮则町=讥心厶）=班勺厶+J + °（工心）（Y ）+0（F ）.ot所以截断误差为：3E=丄 ------ + ---- 「 T h 啰_喟+竺护一 o （F ）T= 0（T + 力”2）解：设点为：（X y ，/林1 ） 3则町=讥勺,_）=以E ，_+1）+ （Y ） +o （巧卩 ot “；：；=班心+1 厶+i ）=叽厶+i ）+滋（ h ）+ * 臥工心）（为 2）+o ox (X)d心;=班心亠心）=班心，/+1）+敕:;D （一力）+ 3 役;D（血 2）+0（亥2）«截断误差为：2舟A 1 ” E= ------------ + ------------ — （―+ _） T h dt dx叭：=班％厶+i ）+敗？心）（_勿+0 @2）〜dx-(史+空八dt dx 呼1_吋】+竺丛Q —O （X ）-(叱 3 +dtdx 22・试用积分插值法推导知铁。
逼近的差分裕式班勺厶叙）一班勺，乩i）+ ——-——£）dtTq2 “-” *\ | (— 4- —)dxdt = | (un t 4- un x)ds = 0* dt & \得-U] /J+U2 r+x^ A-u4 r = 0+JE (j-l? n)F (j,n)G (j^n+l)H (j-l,n+l)^% ~ 的=旳=竹“4 = W/-lMf MTh=h T-T-ll"h + LL r H + ll：4h —LL：N =Op第二章第三章第四章第五章第六章P781.如果①'（0）二0,则称工。
是』（0）的驻点（或稳定:点）-设矩阵A对称（不必正定），求证忑是』（工）的驻点.的充要条件是1心是方程加二&的解B 42・ 试用积分插值法推导知铁。
逼近的差分裕式证： 充分性：①⑻二J 缶)+ 乂(加° -b t ^+—(Ax r x)①'(Ji) = (Ax c - A, x) + A{Ax r x) aEff))S 宀沪若①0)二Q,即(山° 一氛对=0 心怎宀A X Q -h = ()目卩 Ax-b^则帀是方程Ax^b 的解卩 必要性*若心是芳程A^ = b^\解则 Ax a —h - 0 (J 4X 0 — Z?,x) = 0+^◎ (0)=(吐命-b t x) - 0+J所以町是』0)的驻点dpg%3：证明非齐次两点边值间题心現(&)二 e it (E)二 Qu与T 7面的变分间题等价：求血EH 】，认@) = G 使 J(w t ) = min J(y)其中心SiuHU (2)-d』(#) =壬仗站)-(7» —芒⑹戲(D) +而久込叭如(2.13)(提示；先把边值条件齐衩化)+d dxO 字)+梓二/ ax13页证明：令 = w(x) + v(x)其中 w(x) = Q + (x-a)0 w(a) = a yv @) = “v(a) = 0 v(^>) = 0®所以2S = 瞥+qu = j DX DX Pd r /w 血、《, 乂 、 f"丁〔P(T + ：F)]+Q(W + V )" ax dx ax* 丫 d z dv. 产 / d dw 、 豪 令 = - — O —) +(?v = /-(- —^> — +^w) = y;^ ax ax dx ax 所以(1)的等价的形式2厶” =一?0 字)= 卩ax axu(a) = a u\b) = 0a其中久=/-(-£■去字+0W )"ax ax 则由定理22知，讥是辺值间题（2）的解的充要条件是 且满定变分方程"ogf)-C/i 小 0 Vve^Pr (Zv> 一 /j )tdx + p @»： (b)f @) ① W = J(u) = J(u.+^)^— a （u^ + 兔，以.+ 无）一（/,功・ +加）［以・（E ）+加@）］ 2 □2=J(认)+ N[a@・,f)-(/,£)-+乙agd-Qfm 沁卜• Q dx dx 「（加•一/）加x +卩@加：（砂@）-卩@）戊@） Ja（3） => （4）所以可证得• 3必要性：若如 是边值间题（1）的解。
中国海洋大学本科生课程大纲课程属性：公共基础/通识教育/学科基础/专业知识/工作技能，课程性质：必修、选修一、课程介绍1.课程描述：本课程介绍数值求解偏微分方程的基本方法及相关的理论基础。
本课程针对数学类专业高年级（三年级）本科生开设。
课程基本内容包括：有限差分方法、差分格式的稳定性、收敛性分析；变分原理，Galerkin有限元方法等。
通过对模型问题的基本数值方法进行分析，阐明构造数值方法的基本思想和技巧。
通过本课程学习，使学生了解并掌握数值求解偏微分方程的基本思想、基本概念和基本理论（数值格式的相容性、稳定性、收敛性及误差估计等），能够运用算法语言对所学数值方法编制程序在计算机上运行实施并对数值结果进行分析。
培养学生理论联系实际，解决实际问题的能力和兴趣。
2.设计思路：偏微分方程是应用数学的核心内容，在其他科学、技术领域具有广泛深入的应用。
掌握偏微分方程的基础理论及求解方法是数学类专业本科生培养的基本要求。
本课程是在数学物理方程课程基础上开设的延展应用型课程，是一门数值分析理论与实践应用高度融合的专业课。
课程引导学生通过数值方法探讨和理解应用数学工具解决实际- 6 -问题的途径及理论分析框架。
学习本课程需要学生掌握了“数学分析”、“数学物理方程”、“数值分析”及“泛函分析”的核心基本内容。
课程内容安排分为有限差分方法和有限元方法两个单元模块，这是目前应用最广泛、理论发展最完善的两类数值方法，两者既有关联又有本质区别，能够体现偏微分方程数值解法的基本特征。
首先介绍有限差分方法。
有限差分方法是近似求解偏微分方程的应用最广泛的数值方法，以对连续的“导数（微分）”进行离散的“差分”近似为基本出发点，利用Fourier 分析及数值分析的基本理论，讨论椭圆、抛物、双曲等三类典型偏微分方程近似求解方法及近似方法的数学理论分析。
有限元方法是20世纪中期发展起来的基于变分原理的数值方法，具有更直接的物理背景含义，因而受到力学、工程等应用领域广泛的关注和应用。
刘力辉 2010210804011.“画圆为方”问题也是古希腊人所提出几何三大难题中的另一个问题。
即求作一个正方形，使其面积等于已知圆的面积。
不妨设已知圆的半径为 R = 1，试用数值试验显示“画圆为方”问题计算过程中的误差。
（1）MATLAB 程序：y=pi^(1/2); % to generate 15-bit value of square root of pi b=1; d=1; for k=1:8 b=b*10;d=d/10; % b and d combined to control the digit of x x=d*fix(b*y); s(k)=x^3; l(k)=x; endformat long [l',s']（2）误差分析：2． 设，I x xd x n n=+⎰501由 x n = x n + 5 x n – 1 – 5 x n – 1 可得递推式 I n = – 5I n – 1 + 1/ n（1）从 I 0 尽可能精确的近似值出发，利用递推公式：I I nn n =-+-511 ( n = 1,2,…20)计算从 I 1 到 I 20 的近似值；（2）从 I 30 较粗糙的估计值出发，用递推公式：I I nn n -=-+11515 ( n= 30,29, …, 3, 2 )计算从I 1 到 I 20 的近似值；(3) 分析所得结果的可靠性以及出现这种现象的原因。
I 0 =dx x⎰+1051=ln (5+x )10|=ln 6-ln 5 所以I0≈0.18232155679395format longI0=log2(6)/log2(exp(1))-log2(5)/log2(exp(1)) % calculate the value of I0=ln6-ln5 for n=1:20I0=-5*I0+1/n; % recycling equation between I(n+1) and I(n) s(n)=I0; end s'则计算结果为：表1I1 0.0883922160302300 I11 0.0140713362538500 I2 0.0580389198488700 I12 0.0129766520640700 I3 0.0431387340890000 I13 0.0120398166027400 I4 0.0343063295550100 I14 0.0112294884148600 I5 0.0284683522249700 I15 0.0105192245923700 I6 0.0243249055418100 I16 0.0099038770381400 I7 0.0212326151481100 I17 0.0093041442210800 I8 0.0188369242594600 I18 0.0090348344501700 I9 0.0169264898137900 I19 0.0074574066965100 I10 0.0153675509310500I20 0.0127129665174600从计算的数据看出I 20=0.0127129665174600 > I 19=0.0074574066965100又I n 的积分范围为0～1，所以应该有I n >I n+1。
偏微分方程数值方法大作业1、 运用θ-格式编程求解方程：⎪⎪⎪⎩⎪⎪⎪⎨⎧≤≤==≤≤=∂∂≤≤=≤≤≤≤∂∂=∂∂T )t (00t)u(1,t)u(0,1)x (00,(x,0)xu 1)x (0sin πu(x,0)T )t 1,0x (0,x u a tu 22222x解：该双曲型方程为初边值问题，其θ格式为)10(≤≤θ：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+-++==±±===+-++--++--+-+-----+-++-+++-+j j j j j j j k j k j k j kj k j k j k j k j k j k j k j k j a r r a u u j k u u u u u u u u u h a u u u τψϕϕϕϕθθθτ)1()(2,.....2,1........,2,10)]2()2)(21()2([)2(12211221011111111111122112记],,,[121k n k k K u u u U -=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=0110110110 C则A=-2I +C,其中I 为单位矩阵。
由此可把原偏微分方程写成如下矩阵形式：12222222212222])21[(])21()1(4[])21[(-++I -+-+-I -=+I -K K K U C r a r a U C r a r a U C r a r a θθθθθθ化解可得：1]22)21(22)1(4[1]22)2221[(1-+-+-I --+I -=+K U K U C r a r a C r a r a K U θθθθ)1()1(2112112112-⨯-⎥⎥⎥⎥⎥⎥⎦⎥⎢⎢⎢⎢⎢⎢⎣⎢----=N N A取a=1，分别取T=1,2,3，r=2/3,1/4,1/2，h=0.1,0.2,0.3，θ=0,1/4,1/2,1，进行编程计算，根据边界条件可以求得第0层值0j u，再引进虚拟的1-j u，对初始条件进行差分计算，并有t=0时的初始条件满足原偏微分方程，联立可求出第一层的值1u，这样可使边界处的结果也为中心差分，保证了格式的一致精度。
偏微分方程数值解法上机报告（一）一、实验题目：用Ritz-Galerkin 方法求解边值问题2u '',01(0)0,(1)1u x x u u ⎧-+=<<⎨==⎩的第n 次近似()n u x ，基函数()sin(),1,2,...,i x i x i n ϕπ==.二、实验目的：通过本次上机实验，理解求解初值问题的变分问题的最重要的近似解法——Ritz-Galerkin 方法，以便为学习有限元法打好基础。
此外，要熟悉用Matlab 解决数学问题的基本编程方法，提高运用计算机解决问题的能力。
三、实验代码：n=5;syms x;for i=1:np(i)=sin(i*pi*x);q(i)=-i^2*pi^2*sin(i*pi*x);endfor i=1:nb(i)=2*int(p(i),0,1);for j=1:nA(i,j)=int((-q(j)+p(j))*p(i),0,1);endendt=inv(A)*b'四、运行结果：t=2251799813685248/3059521645650671/pi281474976710656/9481460623939047/pi281474976710656/43582901062631895/pi五、总结：通过本次上机，我了解了Ritz-Galerkin 方程 n j j p f cj p i p a n i i ,...,2,1)),(,())(),((1==∑=，明白了用Ritz-Galerkin 方法解决边值问题的变分问题的基本原理，并接近一步提高自己的编程动手能力，受益匪浅。
偏微分方程数值解法上机报告（二）一、 实验题目：用线性元求下列边值问题的数值解2''2sin ,0142(0)0,'(1)0y y x x y y ππ⎧-+=<<⎪⎨⎪==⎩二、 实验目的：通过本次上机，熟悉和掌握用Galerkin 法观点出发导出的求解处置问题数值解的线性有限元法。
偏微分方程的数值解法上机习题四题目：一、 PMECME 算法1) 用四步四阶Adams 方法预估（课本P14，k=3），用三步四阶隐式Adams 方法校正（课本P16，k=3），编写PMECME 算法。
2) 用四步四阶Adams 方法预估（课本P14，k=3），用三步四阶Milne 方法校正（课本P17，（1.2.31）），编写PMECME 算法。
二、 变步计算1) 用误差的后验估计式及三级、四级R-K 方法实现自动选择步长。
计算模型为：()12401u tu u ⎧⎪'=⎨⎪=⎩，（P37，例4.1 ），精确解为()()221u t t =+求解：一、 PMECME 算法1) 用四步四阶Adams 方法预估，用三步四阶隐式Adams 方法校正i. P ： ()1123555937924pn n n n n n hu u f f f f +---=+-+- ii.∵四步四阶Adams 方法的截断误差()()()()()556112511720p n n n u t u h u t O h ++-=+———————————— 三步四阶隐式Adams 方法的截断误差()()()()()55611192720c n n n u t u h u t O h ++-=-+————————————— ∴()()12- 得，()()()55611270720c pn n n u u h u t O h ++-=+ 两边同乘以251720，则有()()()()()556112512513270720c pn n n u u h u t O h ++-=+——两边同乘以19-720，则有()()()()()5561119194270270p cn n n u u h u t O h ++-=-+—— 将()3代入()1得，M1：()11251=+720m pc p n n n n u u u u ++- iii. E1：()111,m mn n n f f t u +++=iv. C ： ()1112919524cm n n n n n n hu u f f f f ++--=++-+ v. 将()4代入()2得，M2：()111119=+720c p cn n n n u u u u ++++- vi.E2：()111,n n n f f t u +++=vii. 编程并作图，将计算解与精确解比较。
ADI隐式交替法三种解法及误差分析（一般的教材上只说第一种）理论部分参看孙志忠：偏微分方程数值解法注意：1.最好不要直接看程序，中间很多公式很烦人的（一定要小心），我写了两天，终于写对了。
2.中间：例如r*(u(i-1,m1,k)+u(i+1,m1,k))形式写成分形式：r*u(i-1,m1,k)+r*u(i+1,m1,k)后面会出错，我也不是很清楚为什么，可能由于舍入误差，或者大数吃掉小数的影响。
3.下面有三个程序4.具体理论看书，先仔细看书（孙志忠：偏微分方程数值解法）或者网上搜一些理论。
Matlab程序：1.function [u u0 p e x y t]=ADI1(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式，截断)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfor i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfor k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))]; u0(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ;endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))]; u0([1 m2+1],j,k)=[exp(0.5*x(j)-t0(k)) exp(0.5*(1+x(j))-t0(k))];endendr=h2/(h1*h1);r1=2*(1-r);r2=2*(1+r);for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组% for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r2*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+...h2*f(i,2,k);for l=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程%for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r2*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);for l=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+... h2*f(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend2.function [u p e x y t]=ADI2(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(D'Yakonov交替方向隐格式)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u'(i,j,k)(引入的过渡层),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;y=(0:m2)*h1+0;t=(0:n)*h2+0;t0=(0:n)*h2+1/2*h2;%定义t0是为了f(x,y,t)~=0的情况% for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));%编程时-t0(k)写成了+t0(k),导致错误;endendend%初始条件for i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endend%边界条件for k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))];endendr=h2/(h1*h1);r4=1+r;r5=r/2;for k=1:nfor i=2:m2u0(i,[1 m1+1],k)=r4*u(i,[1 m1+1],k+1)-r5*(u(i-1,[1 m1+1],... k+1)+u(i+1,[1 m1+1],k+1));endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))];endendr1=r-r*r;r2=2*(r-1)*(r-1);r3=r*r/2;for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组% for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=2*r4*ones(1,m1-1);d(1)=r*u0(i,1,k)+r1*(u(i-1,2,k)+u(i,1,k)+u(i+1,2,k)+...u(i,3,k))+r2*u(i,2,k)+r3*(u(i-1,1,k)+...u(i+1,1,k)+u(i-1,3,k)+u(i+1,3,k))+2*h2*f(i,2,k);for l=2:m1-2d(l)=r1*(u(i-1,l+1,k)+u(i,l,k)+u(i+1,l+1,k)+...u(i,l+2,k))+r2*u(i,l+1,k)+r3*(u(i-1,l,k)+...u(i+1,l,k)+u(i-1,l+2,k)+u(i+1,l+2,k))+2*h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r1*(u(i-1,m1,k)+u(i,m1-1,k)+...u(i+1,m1,k)+u(i,m1+1,k))+r2*u(i,m1,k)+...r3*(u(i-1,m1-1,k)+...u(i+1,m1-1,k)+u(i-1,m1+1,k)+u(i+1,m1+1,k))+2*h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);end%回代过程%u0(i,m1,k)=d(m1-1)/b(m1-1);for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=2*r4*ones(1,m2-1);d(1)=r*u(1,j,k+1)+2*u0(2,j,k);for l=2:m2-2d(l)=2*u0(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=2*u0(m2,j,k)+r*u(m2+1,j,k+1);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend3.function [u u0 p e x y t]=ADI5(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式，未截断)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfor i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfor k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))]; u1(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ;endendr=h2/(h1*h1);r1=2*(1-r);r2=r/4;r3=2*(1+r);for k=1:nfor i=2:m2u0(i,[1 m1+1],k)=u1(i,[1 m1+1],k)-r2*(u(i-1,[1 m1+1],k+1)-...2*u(i,[1 m1+1],k+1)+u(i+1,[1 m1+1],k+1)-u(i-1,[1 m1+1],k)+...2*u(i,[1 m1+1],k)-u(i+1,[1 m1+1],k));endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))];endendfor k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组%for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r3*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+... h2*f(i,2,k);for l=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程%for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r3*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);for l=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+... h2*f(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解 e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend[up e x y t]=ADI2(0.01,0.001,100,100,1000);surf(x,y,e(:,:,1001)) t=1的误差曲面下面是三种方法的误差比较：1.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2）代替）（t=1时的误差）2.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) surf(x,y,e(:,:,11))(表示t=1时的误差）下面是相关数据：1: [u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2)代替）e(:,:,11) =Columns 1 through 60 0 0 0 0 00 0.00040947 0.00025182 0.00019077 0.00017112 0.000176040 0.00057359 0.00042971 0.00035402 0.00032565 0.000336280 0.00066236 0.00054689 0.00047408 0.00044596 0.000462670 0.00072152 0.00062001 0.00055081 0.00052442 0.000545530 0.00076164 0.0006576 0.00058522 0.00055732 0.000579840 0.00078336 0.00065993 0.00057557 0.00054161 0.000562090 0.00078161 0.00061872 0.00051646 0.00047429 0.000489640 0.00073621 0.0005148 0.00039979 0.00035439 0.000363130 0.00056964 0.00031688 0.00022051 0.0001884 0.000191920 0 0 0 0 02.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)e(:,:,11) =Columns 1 through 60 0 0 0 0 00 0.00027006 0.00016305 0.00012104 0.0001071 0.000109950 0.00037754 0.00027817 0.0002253 0.00020483 0.000211160 0.00043539 0.00035386 0.00030207 0.00028124 0.00029140 0.00047398 0.00040104 0.00035113 0.00033111 0.000344050 0.0005003 0.00042535 0.00037309 0.0003519 0.000365710 0.00051479 0.00042699 0.00036681 0.00034164 0.00035410 0.00051415 0.00040056 0.00032887 0.0002985 0.000307640 0.00048504 0.0003335 0.00025411 0.0002221 0.000227060 0.00037609 0.00020532 0.00013956 0.00011718 0.000119020 0 0 0 0 03.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式)e(:,:,11) =Columns 1 through 60 0 0 0 0 00 8.6469e-006 1.4412e-005 1.8364e-005 2.091e-005 2.2174e-0050 1.4412e-005 2.4777e-005 3.2047e-005 3.6716e-005 3.8961e-0050 1.8364e-005 3.2047e-005 4.1789e-005 4.8054e-005 5.1008e-0050 2.091e-005 3.6716e-005 4.8054e-005 5.5353e-005 5.8764e-0050 2.2174e-005 3.8961e-005 5.1008e-005 5.8764e-005 6.2389e-0050 2.2118e-005 3.8698e-005 5.0523e-005 5.8126e-005 6.171e-0050 2.055e-005 3.5581e-005 4.6157e-005 5.2942e-005 5.6197e-0050 1.707e-005 2.8951e-005 3.7128e-005 4.2365e-005 4.4952e-0050 1.0851e-005 1.7698e-005 2.2265e-005 2.5203e-005 2.672e-0050 0 0 0 0 01.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2)代替）Columns 7 through 110 0 0 0 00.00020348 0.00026228 0.00038338 0.00066008 00.00038607 0.00048321 0.00064717 0.00091668 00.00052635 0.00064203 0.00081637 0.0010517 00.0006174 0.00074272 0.00092111 0.0011417 00.00065651 0.00078964 0.00097724 0.0012051 00.00064051 0.00078116 0.00098594 0.0012433 00.00056474 0.00070822 0.00093332 0.0012478 00.00042547 0.00055526 0.00078616 0.0011844 00.00022735 0.00030946 0.00049004 0.00092402 00 0 0 0 02.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)Columns 7 through 110 0 0 0 00.00012826 0.00016798 0.00024986 0.00043637 00.00024444 0.00031023 0.00042173 0.00060513 00.00033401 0.00041257 0.00053179 0.00069358 00.00039216 0.00047742 0.00059986 0.00075263 00.00041704 0.00050761 0.00063642 0.00079439 00.00040657 0.0005021 0.00064226 0.00081984 00.00035784 0.000455 0.00060828 0.00082334 00.00026866 0.00035628 0.00051263 0.0007824 00.00014262 0.00019789 0.00031956 0.0006113 00 0 0 0 0 3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) Columns 7 through 110 0 0 0 02.2118e-005 2.055e-005 1.707e-005 1.0851e-005 03.8698e-005 3.5581e-005 2.8951e-005 1.7698e-005 05.0523e-005 4.6157e-005 3.7128e-005 2.2265e-005 05.8126e-005 5.2942e-005 4.2365e-005 2.5203e-005 06.171e-005 5.6197e-005 4.4952e-005 2.672e-005 06.1116e-005 5.5803e-005 4.4817e-005 2.6785e-005 05.5803e-005 5.1239e-005 4.1529e-005 2.5153e-005 04.4817e-005 4.1529e-005 3.42e-005 2.126e-005 02.6785e-005 2.5153e-005 2.126e-005 1.3869e-005 00 0 0 0 0。
数值计算方法试题一一、 填空题（每空1分，共17分）1、如果用二分法求方程043=-+x x 在区间]2,1[内的根精确到三位小数，需对分（ ）次。
2、迭代格式)2(21-+=+k k k x x x α局部收敛的充分条件是α取值在（ ）。
3、已知⎪⎩⎪⎨⎧≤≤+-+-+-≤≤=31)1()1()1(2110)(233x c x b x a x x x x S 是三次样条函数，则a =( )，b =（ ），c =（ ）。
4、)(,),(),(10x l x l x l n 是以整数点n x x x ,,,10 为节点的Lagrange 插值基函数，则∑==nk kx l0)(( )，∑==nk k jk x lx 0)(( )，当2≥n 时=++∑=)()3(204x l x xk k n k k( )。
5、设1326)(247+++=x x x x f 和节点,,2,1,0,2/ ==k k x k 则=],,,[10n x x x f 和=∆07f。
6、5个节点的牛顿-柯特斯求积公式的代数精度为 ，5个节点的求积公式最高代数精度为 。
7、{}∞=0)(k kx ϕ是区间]1,0[上权函数x x =)(ρ的最高项系数为1的正交多项式族，其中1)(0=x ϕ，则⎰=14)(dx x x ϕ 。
8、给定方程组⎩⎨⎧=+-=-221121b x ax b ax x ，a 为实数，当a 满足 ，且20<<ω时，SOR 迭代法收敛。
9、解初值问题00(,)()y f x y y x y '=⎧⎨=⎩的改进欧拉法⎪⎩⎪⎨⎧++=+=++++)],(),([2),(]0[111]0[1n n n n n n n n n n y x f y x f h y y y x hf y y 是阶方法。
10、设⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=11001a a a a A ，当∈a （ ）时，必有分解式T LL A =，其中L 为下三角阵，当其对角线元素)3,2,1(=i l ii 满足（ ）条件时，这种分解是唯一的。
一维抛物线偏微分方程数值解法（2）上一篇文章请参看一维抛物线偏微分方程数值解法（1）解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法）Ut-Uxx=0, 0<x<1,0<t<=1（Ut-aUxx=f(x,t),a>0)U(x,0)=e^x, 0<=x<=1,U(0,t)=e^t,U(1,t)=e^(1+t), 0<t<=1精确解为：U(x,t)=e^(x+t);Matlab程序：（此为向后差分法）function [u p e x t]=pwxywxh(h1,h2,m,n)%欧拉向后差分法解一维抛物线型偏微分方程%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m,n分别为空间，时间网格数%p为精确解，u为数值解,e为误差x=(0:m)*h1+0;t=(0:n)*h2+0;for(i=1:n+1)for(j=1:m+1)f(i,j)=0;endendfor(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endr=h2/(h1*h1);for(i=2:n+1) %外循环，先固定每一时间层，每一时间层上解一线性方程组%a(1)=0;b(1)=1+2*r;c(1)=-r;d(1)=u(i-1,2)+h2*f(i,2)+r*u(i,1);for(k=2:m-2)a(k)=-r;b(k)=1+2*r;c(k)=-r;d(k)=u(i-1,k+1)+h2*f(i,k+1);%输入部分系数矩阵，为0的矩阵元素不输入%enda(m-1)=-r;b(m-1)=1+2*r;d(m-1)=u(i-1,m)+h2*f(i,m)+r*u(i,m+1);for(k=1:m-2) %开始解线性方程组消元过程a(k+1)=-a(k+1)/b(k);b(k+1)=b(k+1)+a(k+1)*c(k);d(k+1)=d(k+1)+a(k+1)*d(k);endu(i,m)=d(m-1)/b(m-1); %回代过程%for(k=m-2:-1:1)u(i,k+1)=(d(k)-c(k)*u(i,k+2))/b(k);endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j)+t(i)); %p为精确解e(i,j)=abs(u(i,j)-p(i,j));%e为误差endend[u p e x t]=pwxywxh(0.1,0.005,10,200);surf(x,t,e);xlabel('x');ylabel('t');zlabel('e');>> title('误差曲面');plot(t,e)误差较之前的欧拉向前差分格式增长了两倍[u p e x t]=pwxywxh(0.1,0.05,10,20); plot(t,e)[u p e x t]=pwxywxh(0.01,0.05,100,20); plot(t,e)[u p e x t]=pwxywxh(0.01,0.005,100,200);plot(x,e)[u p e x t]=pwxywxh(0.005,0.005,200,200); plot(x,e)X=1时，出现了误差？？? 不是边界条件吗？不能理解这方法还是比前一种方法误差大呀不过可以随便改变时间、空间步长。
偏微分方程数值解——第一次大作业作业一：用有限差分法计算：()()()[]()()-,0,00d du a x c x u f x x L dx dx u u L ⎧⎛⎫+=∈⎪ ⎪⎝⎭⎨⎪==⎩要求：1）L 随机选择2）()(),a x c x 光滑（随机），且()()min max 0,0a a x a c x <≤≤<∞≥ 3）统计L 2-误差 4）规模尽可能大 根据题意：将线段0~L 均匀分成M 段，则网格步长/h L M =，节点分布为0,1,2,,M 共M+1个，k x kh =，可得原方程的差分形式：1111220,1,2,,10k k k k k k k k k M u u u u a a h h c u f k M h u u +-+---⎧-⎪⎪-+==-⎨⎪==⎪⎩ 其中，下标k 表示在k x 处取值，12k +（12k -）表示在k x 与1k x +（k x 与1k x -）中点处的取值。
化简上式可得：1111222211222,1,2,,1k k k k k k k k k aa a a u c u u f k M h h h --++-++⎛⎫⎪-++-==- ⎪ ⎪⎝⎭写成矩阵形式：AU g =，得()()13322212233552222222252253322222233122212211000000M M M M M M M M M M a a a c h h a a a a c h hha A haa a c hha a a c h h ---------⨯-+⎛⎫ ⎪+-⎪⎪+ ⎪ ⎪-+- ⎪ ⎪ ⎪=-⎪⎪+ ⎪ ⎪+- ⎪ ⎪+ ⎪ ⎪-+⎝⎭12211M M M u u U u u ---⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭ ，12211M M M f f g f f ---⎛⎫⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭计算时，取1L =，()x a x e =，()2x c x e =，()2f x =，可以解得原方程的解析解（精确解）为()3233111x x xe e e u x e e e e e ----=++--由于A 为对称三对角阵且主对角占优，可以用直接法求解此方程组。
有限差分法求解偏微分方程摘要：本文主要使用有限差分法求解计算力学中的系统数学模型，推导了有限差分法的理论基础，并在此基础上给出了部分有限差分法求解偏微分方程的算例验证了推导的正确性及操作可行性。
关键词：计算力学，偏微分方程，有限差分法Abstract：This dissertation mainly focuses on solving the mathematic model of computation mechanics with finite-difference method. The theoretical basis of finite-difference is derived in the second part of the dissertation, and then I use MATLAB to program the algorithms to solve some partial differential equations to confirm the correctness of the derivation and the feasibility of the method.Key words：Computation Mechanics, Partial Differential Equations, Finite-Difference Method1 引言机械系统设计常常需要从力学观点进行结构设计以及结构分析，而这些分析的前提就是建立工程问题的数学模型。
通过对机械系统应用自然的基本定律和原理得到带有相关边界条件和初始条件的微分积分方程，这些微分积分方程构成了系统的数学模型。
求解这些数学模型的方法大致分为解析法和数值法两种，而解析法的局限性众所周知，当系统的边界条件和受载情况复杂一点，往往求不出问题的解析解或近似解。
另一方面，计算机技术的发展使得计算更精确、更迅速。
偏微风方程数值解法第2章 有限差分方法的基本概念习题(p.44)第二题讨论对流方程0u ut x∂∂+=∂∂,a>0 的差风格式 11110n n n n j jj ju u u u ahτ++++--+=的截断误差及稳定解：(1)先计算截断误差：用Taylor 级数在n ju 展开： (,)(,)(,)(,)11()()x t t x t t x x t t x t T u a u u h∆∆∆∆τ++++=+-＝(,)(,)(,)(,)(,)(,)1()()x t t x t x x t t x t x t t x t a u u u u u u h∆∆∆∆τ++++-+--+=2322332322232321111()26262u u a u u u ah uh a ah t t h x x t x x t ττττ∂∂∂∂∂∂++--++∂∂∂∂∂∂∂∂= +∂∂+∂∂222221xuah t u τ 可见：),(0),(),(),(τττ+∂∂=-+tt x u t x u t x u(,)(,)(,)0()u x h t u x t u x t h h xττ++-+∂=+∂所以，截断误差为()o h τ+ （2）稳定性：1111()n n n n j j j j u u a u u λ++++=+-将ikjhn n j ev u =代入上式得： 11(1)1()n ikjh n ikjh n ik j h n ikhj v e v e a v e v e λ++++=+- n v =))1(1(1ikh n e a v -+-+λ 得增长因子；1(,)1(1)ikh G k a e τλ=+-=11(cos 1)sin a kh ia khλλ+-+[]222221(,)1(cos 1)sin G k a kh a khτλλ=+-+=22111222(1)cos a a a a khλλλλ≤+-+-0>a 时对于任何值2(,)1k G τ≤，所以此格式无条件稳定。
偏微分方程数值解大作业目录第一题 (3)第二题 (7)第三题 (16)第四题 (20)第五题 (26)第六题（附加题1） (39)第七题（附加题2） (45)第八题（附加题3） (51)第一题习题13.（1）解曲线图图1 （2）误差曲线图图2（3）表格表1 部分点处精确解和取不同步长时所得的数值解表2 取不同步长时部分结点处数值解的误差的绝对值和数值解的最大误差（4）MATLAB源代码M=64;a=0;b=pi/2;h=(b-a)/M;x=[a+h:h:b-h];u=zeros(M-1,M-1);u(1,1)=(2/h^2)+(x(1)-1/2)^2;u(1,2)=-(1/h^2);u(M-1,M-1)=(2/h^2)+(x(M-1)-1/2)^2;u(M-1,M-2)=-(1/h^2);for i=2:M-2u(i,i-1)=-(1/h^2);u(i,i)=(2/h^2)+(x(i)-1/2)^2;u(i,i+1)=-(1/h^2);endf=zeros(M-1,1)f(1)=(x(1).*x(1)-x(1)+5/4).*sin(x(1));f(M-1)=(x(M-1).*x(M-1)-x(M-1)+5/4).*sin(x(M-1))+1/h^2; for j=2:M-2f(j)=(x(j).*x(j)-x(j)+5/4).*sin(x(j));endy=inv(u)*f; true=sin(x); plot(x,y'-true)第二题习题二（P67 第3题）（1）h=1/4, τ=1/4精确解数值解误差（2）h=1/8, τ=1/8精确解数值解误差（3）h=1/16, τ=1/16精确解数值解误差（4）h=1/32, τ=1/32精确解数值解误差（5）h=1/64, τ=1/64精确解精确解误差（6）表格（7）Matlab代码function[p,e,u,x,y,k]=fivepoint(h,m,n,kmax,ep)%五点差分法和G-S迭代法解椭圆型方程%kmax为最大迭代次数；%m,n分别为x,y方向的网格数；%ep为精度；%u为差分解，p为精确解，e为误差；%例如在命令行窗口输入[p,e,u,x,y,k]=fivepoint(1/64,64,64,10000,1e-10);%代表步长为1/64,精度为10^(-10),最大迭代次数为10000的五点差分格式syms temp;u=zeros(n+1,m+1);x=0+(0:m)*h;y=0+(0:n)*h;for(j=1:n+1)u(j,1)=sin(y(j))+cos(y(j));u(j,m+1)=exp(1)*(sin(y(j))+cos(y(j)));endfor(i=1:m+1)u(1,i)=exp(x(i));u(n+1,i)=exp(x(i))*(sin(1)+cos(1));endt=zeros(n-1,m-1);for(k=1:kmax)for(j=2:n)for(i=2:m)temp=(u(j,i+1)+u(j,i-1)+u(j+1,i)+u(j-1,i))/4; t(j,i)=(temp-u(j,i))*(temp-u(j,i));u(j,i)=temp;endendt(j,i)=sqrt(t(j,i));if(k>kmax)break;endif(max(max(t))<ep)break;endendfor(j=1:n+1)for(i=1:m+1)p(j,i)=(sin(y(j))+cos(y(j)))*exp(x(i));e(j,i)=abs(u(j,i)-p(j,i));endend代码使用说明：在命令行窗口输入[p,e,u,x,y,k]=fivepoint(1/64,64,64,10000,1e-10);代表步长为1/64,精度为10^(-10),最大迭代次数为10000的五点差分格式surf(x,y,p)可作出近似解曲线surf(x,y,u)可作出精确解曲线surf(x,y,e)可作出误差曲线注意精度不要选取的太低，否则随着步长减小误差反而增大。
三张图分别为步长为1/4,1/16,1/64的误差曲线（精度均为10^(-10)）。
第三题习题三P138 第12题（1）h=0.010 tau=0.010k (x,t) 数值解精确解误差10 (0.5,0.1) 0.641417 0.642042 6.256223e-0420 (0.5,0.2) 0.486439 0.487230 7.914001e-0430 (0.5,0.3) 0.326667 0.327550 8.834544e-0440 (0.5,0.4) 0.163640 0.164597 9.575961e-0450 (0.5,0.5) -0.001021 0.000000 1.020916e-0360 (0.5,0.6) -0.165671 -0.164597 1.073862e-0370 (0.5,0.7) -0.328666 -0.327550 1.116054e-0380 (0.5,0.8) -0.488378 -0.487230 1.147092e-0390 (0.5,0.9) -0.643209 -0.642042 1.166668e-03100 (0.5,1.0) -0.791614 -0.790439 1.174587e-03h=0.005 tau=0.005k (x,t) 数值解精确解误差20 (0.5,0.1) 0.641729 0.642042 3.129321e-0440 (0.5,0.2) 0.486834 0.487230 3.959919e-0460 (0.5,0.3) 0.327108 0.327550 4.420366e-0480 (0.5,0.4) 0.164118 0.164597 4.790797e-04100 (0.5,0.5) -0.000511 0.000000 5.107000e-04120 (0.5,0.6) -0.165135 -0.164597 5.371294e-04140 (0.5,0.7) -0.328109 -0.327550 5.581797e-04160 (0.5,0.8) -0.487804 -0.487230 5.736512e-04180 (0.5,0.9) -0.642626 -0.642042 5.833907e-04200 (0.5,1.0) -0.791026 -0.790439 5.873012e-04h=0.0050 tau=0.0005k (x,t) 数值解精确解误差200 (0.5,0.1) 0.642011 0.642042 3.115887e-05400 (0.5,0.2) 0.487191 0.487230 3.948351e-05600 (0.5,0.3) 0.327506 0.327550 4.412446e-05800 (0.5,0.4) 0.164550 0.164597 4.786762e-051000 (0.5,0.5) -0.000051 0.000000 5.106902e-051200 (0.5,0.6) -0.164651 -0.164597 5.375135e-05 1400 (0.5,0.7) -0.327606 -0.327550 5.589538e-05 1600 (0.5,0.8) -0.487288 -0.487230 5.748076e-05 1800 (0.5,0.9) -0.642101 -0.642042 5.849178e-05 2000 (0.5,1.0) -0.790498 -0.790439 5.891837e-05h τ E∞(h, τ) E∞(2h,2τ)/E∞(h, τ)1/10 1/10 1.175210e-02 *1/20 1/20 5.910745e-03 1.9881/40 1/40 2.955702e-03 1.9991/80 1/80 1.478257e-03 1.9991/160 1/160 7.391688e-04 1.999 (2) t=1时数值解的误差曲线(3)外推步长1/100 步长1/200 数值解精确解误差（3）Matlab代码：h=input('输入空间步长h=');tau=input('输入时间步长tau=');%h=1/100;tau=1/200;r=tau/h^2;m=1/h; n=1/tau;i=1:m-1; k=1:n-1;A=diag(-r*ones(m-2,1),-1)+diag((1+2*r)*ones(m-1,1))+diag(-r*ones(m-2,1),1);B=diag(r*ones(m-2,1),-1)+diag((1-2*r)*ones(m-1,1))+diag(r*ones(m-2,1),1);%C=diag(-1/2*s^2*ones(m-2,1),-1)+diag((2+s^2)*ones(m-1,1))+diag(-1/2*s^2*ones(m-2,1),1);xi=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;ui0=(exp(xi)*sin(1/2))'; u00=sin(1/2); um0=exp(1)*sin(1/2);u0k=sin(1/2-tk); umk=exp(1)*sin(1/2-tk);fik=zeros(m-1,n);for j=1:nfik(:,j)=-exp(xi)*(cos(1/2-tk(j))+2*sin(1/2-tk(j)));enduik=zeros(m-1,n);uik(:,1)=inv(A)*(B*ui0+[r*u00;zeros(m-3,1);r*um0]+[r*u0k(1);zeros(m-3,1);r*umk(1)]+tau*(-exp(xi)*(cos(1/2)+2*sin(1/2)))');%first k=1%uik(:,2)=inv(A)*(B*uik(:,1)+[r*u0k(1);zeros(m-3,1);r*umk(1)]+[r*u0k(2);zeros(m-3,1);r*umk(2)]+tau*fik(:,1));%%%%%%%%%%%%%%%%for k=2:n-1D=(B*uik(:,k)+[r*u0k(k);zeros(m-3,1);r*umk(k)]+[r*u0k(k+1);zeros(m-3,1);r*umk(k+1)]+tau*fik(:,k));uik(:,k+1)=inv(A)*D;endtrue=zeros(m-1,n);for i=1:m-1true(i,:)=exp(xi(i))*sin(1/2-tk);endfprintf('h=%.4f tau=%.4f\n',h,tau);fprintf('k (x,t) 数值解精确解误差\n');for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));elsefprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));endenda=max(abs(uik-true));E=max(a);fprintf('Emax=%e',E);[x,t]=meshgrid(xi,tk);mesh(x,t,abs(uik-true)');xlabel('xi');ylabel('tk');hold on;%e=max(true-uik);% h=max(e);% plot(e)%surf(abs(uik-true))%surf(true)% plot(xi,abs(uik(:,n)-true(:,n)))% hold on;第四题P176 习题四第8题(1)表格h=0.010 tau=0.005(x,t) 数值解精确解误差(0.5,0.1) 0.049979 0.049979 4.791843e-08 (0.5,0.2) 0.099833 0.099833 7.083966e-08 (0.5,0.3) 0.149438 0.149438 4.406409e-08 (0.5,0.4) 0.198669 0.198669 5.641112e-08 (0.5,0.5) 0.247404 0.247404 2.789924e-07 (0.5,0.6) 0.295519 0.295520 9.784169e-07 (0.5,0.7) 0.342896 0.342898 1.775191e-06 (0.5,0.8) 0.389416 0.389418 2.623113e-06 (0.5,0.9) 0.434962 0.434966 3.474007e-06 (0.5,1.0) 0.479421 0.479426 4.279944e-06h=0.005 tau=0.010(x,t) 数值解精确解误差(0.5,0.1) 0.049979 0.049979 1.917987e-07(0.5,0.2) 0.099834 0.099833 2.836035e-07 (0.5,0.3) 0.149438 0.149438 1.765491e-07 (0.5,0.4) 0.198669 0.198669 2.261917e-07 (0.5,0.5) 0.247403 0.247404 1.148338e-06 (0.5,0.6) 0.295516 0.295520 3.962614e-06 (0.5,0.7) 0.342891 0.342898 7.167473e-06 (0.5,0.8) 0.389408 0.389418 1.059995e-05 (0.5,0.9) 0.434951 0.434966 1.405730e-05 (0.5,1.0) 0.479408 0.479426 1.748442e-05h τE∞(h, τ) E∞(2h,2τ)/E∞(h, τ)1/10 1/10 1.767294e-03 *1/20 1/20 4.356573e-04 4.057 1/40 1/40 1.100309e-04 3.959 1/80 1/80 2.753878e-05 3.995 1/160 1/160 6.882379e-06 4.001 1/320 1/320 1.720007e-06 4.001(2)精确解 (h=1/100, τ=1/100)数值解(h=1/10, τ=1/10)(3)上曲面 (h=1/10, τ=1/10) 中曲面(h=1/20, τ=1/20) 下曲面(h=1/40, τ=1/40)（4）Matlab代码：h=input('输入空间步长h=');tau=input('输入时间步长tau='); s=tau/h;m=1/h; n=1/tau;i=1:m-1; k=1:n-1;A=diag(-1/2*s^2*ones(m-2,1),-1)+diag((1+s^2)*ones(m-1,1))+diag(-1/2*s^2*ones(m-2,1),1);B=diag(1/2*s^2*ones(m-2,1),-1)+diag(-(1+s^2)*ones(m-1,1))+diag(1/2*s^2*ones(m-2,1),1);C=diag(-1/2*s^2*ones(m-2,1),-1)+diag((2+s^2)*ones(m-1,1))+diag(-1/2*s^2*ones(m-2,1),1);xi=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;ui0=zeros(m-1,1); u00=0; um0=0;u0k=zeros(n,1); umk=sin(tk);fik=zeros(m-1,n);for j=1:nfik(:,j)=(((tk(j)).^2-(xi).^2).*sin((xi).*(tk(j))))'; enduik=zeros(m-1,n);uik(:,1)=inv(C)*(2*tau*xi'+[zeros(m-2,1);1/2*s^2*umk(1)]);%first k=1%D=2*uik(:,1)+B*ui0+tau^2*fik(:,1)+[zeros(m-2,1);1/2*s^2*umk(2)];uik(:,2)=inv(A)*D;%%%%%%%%%%%%%%%%for k=2:n-1D=2*uik(:,k)+B*uik(:,k-1)+[zeros(m-2,1);1/2*s^2*umk(k-1)]+tau^2*fik(:,k)+[zeros(m-2,1);1/2*s^2*umk(k+1)];uik(:,k+1)=inv(A)*D;endtrue=zeros(m-1,n);for i=1:m-1true(i,:)=sin(xi(i)*tk);endfprintf('h=%.3f tau=%.3f\n',h,tau);fprintf('k (x,t) 数值解精确解误差\n');for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),tru e(m/2,p),abs(uik(m/2,p)-true(m/2,p)));elsefprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),tru e(m/2,p),abs(uik(m/2,p)-true(m/2,p)));endenda=max(abs(uik-true));E=max(a);fprintf('Emax=%e',E);[x,t]=meshgrid(xi,tk);mesh(x,t,true');xlabel('xi');ylabel('tk');hold on;%e=max(true-uik);% h=max(e)% plot(e)%surf(abs(uik-true))%surf(true)%plot(abs(uik(:,n)-true(:,n)))第五题习题56.（1）曲线图① h=0.05 τ= 0.05数值解精确解误差曲线图② h=0.025 τ= 0.025数值解精确解误差曲线图（2）表格① h=0.0500 τ=0.0500k (x,y,t) 数值解精确解误差02 (0.5,0.5,0.1) 0.022374 0.024997 2.623812e-03 04 (0.5,0.5,0.2) 0.049668 0.049979 3.116082e-04 06 (0.5,0.5,0.3) 0.074901 0.074930 2.839556e-05 08 (0.5,0.5,0.4) 0.099852 0.099833 1.849222e-05 10 (0.5,0.5,0.5) 0.124713 0.124675 3.871957e-05 12 (0.5,0.5,0.6) 0.149497 0.149438 5.852752e-05 14 (0.5,0.5,0.7) 0.174189 0.174108 8.110513e-05 16 (0.5,0.5,0.8) 0.198776 0.198669 1.066039e-04 18 (0.5,0.5,0.9) 0.223241 0.223106 1.347368e-04 20 (0.5,0.5,1.0) 0.247569 0.247404 1.651274e-04② h=0.0250 τ=0.0250k (x,y,t) 数值解精确解误差04 (0.5,0.5,0.1) 0.024101 0.024997 8.967699e-04 08 (0.5,0.5,0.2) 0.049855 0.049979 1.240867e-04 12 (0.5,0.5,0.3) 0.074916 0.074930 1.412105e-05 16 (0.5,0.5,0.4) 0.099837 0.099833 3.688824e-06 20 (0.5,0.5,0.5) 0.124684 0.124675 9.652543e-06 24 (0.5,0.5,0.6) 0.149453 0.149438 1.475416e-05 28 (0.5,0.5,0.7) 0.174129 0.174108 2.044073e-05 32 (0.5,0.5,0.8) 0.198696 0.198669 2.684217e-05 36 (0.5,0.5,0.9) 0.223140 0.223106 3.389905e-05 40 (0.5,0.5,1.0) 0.247445 0.247404 4.151868e-05（3）MATLAB源代码clear;%%h=input('输入空间步长h=');tau=input('输入时间步长tau=');r=tau/h^2;m=1/h;n=1/tau;i=1:m-1;j=1:m-1;k=1:n-1;A=diag(-r/2*ones(m-2,1),-1)+diag((1+r)*ones(m-1,1))+diag(-r/2*ones(m-2,1),1);B=diag(r/2*ones(m-2,1),-1)+diag((1-r)*ones(m-1,1))+diag(r/2*ones(m-2,1),1); xi=0+1/m:1/m:1-1/m;yj=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;uij0=zeros(m-1,m-1);u000=0;u0m0=0;um00=0;umm0=0;uim0=zeros(m-1,1);um0k=0;ummk=sin(tk);umj0=0;ui00=zeros(m-1,1);%%u0jk=zeros(m-1,n);umjk=zeros(m-1,n);ui0k=zeros(m-1,n);uimk=zeros(m-1,n);for k=1:nu0jk(:,k)=0;umjk(:,k)=sin(yj*tk(k))';ui0k(:,k)=0;uimk(:,k)=sin(xi*tk(k))';endfijk=zeros(m-1,m-1,n);fij0=0;for k=1:nfor j=1:m-1fijk(:,j,k)=((xi.^2+yj(j).^2).*(tk(k).^2).*sin(xi*yj(j)*tk(k))+xi.*yj(j).*cos( xi*yj(j)*tk(k)))';endenduijk0p5=zeros(m-1,m-1,n);uijk=zeros(m-1,m-1,n);true=zeros(m-1,m-1,n);for k=1:nfor j=1:m-1true(:,j,k)=sin(xi*yj(j)*tk(k))';endend%%for k=1:nfor j=1:m-1if j==1if k==1u0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*0))-1/4*tau*1/h^2*((um0k-2*umjk(j,k)+umjk(j+1,k)-(um0k-2*umj0+umj0)));elseu0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((um0k-2*umjk(j,k)+umjk(j+1,k)-(um0k-2*umjk(j,k-1)+umjk(j+1,k-1))));endelseif j==m-1if k==1u0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+ummk(k)-(umj0-2*umj0+umm0)));elseu0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+ummk(k)-(umjk(j-1,k-1)-2*umjk(j,k-1)+ummk(k-1))));endelseif k==1u0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+umjk(j+1,k)-(umj0-2*umj0+umj0)));elseu0jk0p5=0;umjk0p5=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+umjk(j+1,k)-(umjk(j-1,k-1)-2*umjk(j,k-1)+umjk(j+1,k-1))));endendif k==1if j==1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uij0(p,1);if p~=1temp1(p,p-1)=uij0(p-1,2);endendtemp2=ui00+[zeros(m-2,1);uij0(m-1,2)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fij0)/2;temp4=diag((1-r)*ones(m-1,1))+diag(r/2*ones(m-2,1),1);uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); elseif j==m-1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uij0(p,m-2);if p~=1temp1(p,p-1)=uij0(p-1,m-1);endendtemp2=uim0+[zeros(m-2,1);2*(1-r)/r*uij0(m-1,m-1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fij0)/2;temp4=diag(r/2*ones(m-1,1))+diag((1-r)*ones(m-2,1),1);uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); elsetemp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uij0(p,j);if p~=1temp1(p,p-1)=uij0(p-1,j+1);endif p~=m-1temp1(p,p+1)=uij0(p+1,j-1);endendtemp2=[uij0(1,j-1);zeros(m-3,1);uij0(m-1,j+1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fij0)/2;temp4=B;uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); endelseif j==1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk(p,1,k-1);if p~=1temp1(p,p-1)=uijk(p-1,2,k-1);endendtemp2=ui0k(:,k-1)+[zeros(m-2,1);uijk(m-1,2,k-1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fijk(:,j,k-1))/2;temp4=diag((1-r)*ones(m-1,1))+diag(r/2*ones(m-2,1),1);uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3);elseif j==m-1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk(p,m-2,k-1);if p~=1temp1(p,p-1)=uijk(p-1,m-1,k-1);endendtemp2=uimk(:,k-1)+[zeros(m-2,1);2*(1-r)/r*uijk(m-1,m-1,k-1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%all add in temp2temp3=(fijk(:,j,k)+fijk(:,j,k-1))/2;temp4=diag(r/2*ones(m-1,1))+diag((1-r)*ones(m-2,1),1);uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); elsetemp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk(p,j,k-1);if p~=1temp1(p,p-1)=uijk(p-1,j+1,k-1);endif p~=m-1temp1(p,p+1)=uijk(p+1,j-1,k-1);endendtemp2=[uijk(1,j-1,k-1);zeros(m-3,1);uijk(m-1,j+1,k-1)]+[u0jk0p5;zeros(m-3,1);umjk0p5];%%%%%%%%%%%%%%%%%%addtemp3=(fijk(:,j,k)+fijk(:,j,k-1))/2;temp4=B;uijk0p5(:,j,k)=A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3); endendend%%%%%%%%%%%%%%%%%%%%%%edge Vec%%%%%%%%%%%%%%%%%%%u0jk0p5vec=0;umjk0p5vec=zeros(m-1,n);%%used before and now is freeif k==1for j=1:m-1if j==1u0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((um0k-2*umjk(j,k)+umjk(j+1,k)-(umj0-2*umj0+umj0)));elseif j==m-1u0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+ummk(k)-(umj0-2*umj0+umj0)));elseu0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+0)-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+umjk(j+1,k)-(umj0-2*umj0+umj0))); endendelsefor j=1:m-1if j==1u0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((um0k-2*umjk(j,k)+umjk(j+1,k)-(um0k-2*umjk(j,k-1)+umjk(j+1,k-1))));elseif j==m-1u0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+ummk(k)-(umjk(j-1,k-1)-2*umjk(j,k-1)+ummk(k-1))));elseu0jk0p5vec=0;umjk0p5vec(j,k)=1/2*(sin(yj(j)*tk(k))+sin(yj(j)*tk(k-1)))-1/4*tau*1/h^2*((umjk(j-1,k)-2*umjk(j,k)+umjk(j+1,k)-(umjk(j-1,k-1)-2*umjk(j,k-1)+umjk(j+1,k-1))));endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for i=1:m-1if i==1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk0p5(1,p,k);if p~=1temp1(p,p-1)=uijk0p5(2,p-1,k);endendtemp2=u0jk0p5vec+[zeros(m-2,1);uijk0p5(2,m-1,k)]+[ui0k(1,k);zeros(m-3,1);uimk(1,k)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if k==1temp3=(fijk(i,:,k)+fij0)'/2;elsetemp3=(fijk(i,:,k)+fijk(i,:,k-1))'/2;endtemp4=diag((1-r)*ones(m-1,1))+diag(r/2*ones(m-2,1),1);uijk(i,:,k)=(A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3))';elseif i==m-1temp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk0p5(m-2,p,k);if p~=1temp1(p,p-1)=uijk0p5(m-1,p-1,k);endendtemp2=umjk0p5vec(:,k)+[zeros(m-2,1);(1-r)*2/r*uijk0p5(m-1,m-1,k)]+[ui0k(m-1,k);zeros(m-3,1);uimk(m-1,k)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if k==1temp3=(fijk(i,:,k)+fij0)'/2;elsetemp3=(fijk(i,:,k)+fijk(i,:,k-1))'/2;endtemp4=diag((r/2)*ones(m-1,1))+diag((1-r)*ones(m-2,1),1);uijk(i,:,k)=(A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3))';elsetemp1=zeros(m-1,m-1);for p=1:m-1temp1(p,p)=uijk0p5(i,p,k);if p~=1temp1(p,p-1)=uijk0p5(i+1,p-1,k);endif p~=m-1temp1(p,p+1)=uijk0p5(i-1,p+1,k);endendtemp2=[uijk0p5(i-1,1,k);zeros(m-3,1);uijk0p5(i+1,m-1,k)]+[ui0k(i,k);zeros(m-3,1);uimk(i,k)];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if k==1temp3=(fijk(i,:,k)+fij0)'/2;elsetemp3=(fijk(i,:,k)+fijk(i,:,k-1))'/2;endtemp4=B;uijk(i,:,k)=(A\(diag(temp4*temp1)+r/2*temp2+tau/2*temp3))';endendendsurf(uijk(:,:,n));figure(2)surf(true(:,:,n));figure(3)surf(uijk(:,:,n)-true(:,:,n));fprintf('h=%.4f tau=%.4f\n',h,tau);fprintf('k (x,y,t) 数值解精确解误差 \n');for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,0.5,%.1f) %f %f %e\n',p,0.1*k,uijk(m/2,m/2,p),true(m/2,m/2 ,p),abs(uijk(m/2,m/2,p)-true(m/2,m/2,p)));elsefprintf('%d(0.5,0.5,%.1f) %f %f %e\n',p,0.1*k,uijk(m/2,m/2,p),true(m/2,m/2 ,p),abs(uijk(m/2,m/2,p)-true(m/2,m/2,p)));endend第六题补充题1（教材P94/算例3.4）（1）差分格式ðu ðt (x i,t k)−að2uðt2(x i,t k)=f(x i,t k)u i k+1+u i k−12τ−aδx2u i k−1+δx2u i k+12=f(x i,t k)（2）算例3.4部分结点数值解精确解和误差的绝对值① h=0.100 τ=0.010k (x,t) 数值解精确解误差10 (0.5,0.1) 1.822237 1.822119 1.181012e-0420 (0.5,0.2) 2.013930 2.013753 1.773622e-0430 (0.5,0.3) 2.225754 2.225541 2.135582e-0440 (0.5,0.4) 2.459846 2.459603 2.425890e-0450 (0.5,0.5) 2.718552 2.718282 2.705635e-0460 (0.5,0.6) 3.004466 3.004166 2.999407e-0470 (0.5,0.7) 3.320449 3.320117 3.318310e-0480 (0.5,0.8) 3.669664 3.669297 3.668593e-0490 (0.5,0.9) 4.055605 4.055200 4.054907e-04100 (0.5,1.0) 4.482137 4.481689 4.481546e-04② h=0.050 τ=1/400k (x,t) 数值解精确解误差40 (0.5,0.1) 1.822148 1.822119 2.874976e-05 80 (0.5,0.2) 2.013796 2.013753 4.313840e-05 120 (0.5,0.3) 2.225593 2.225541 5.191899e-05 160 (0.5,0.4) 2.459662 2.459603 5.896391e-05 200 (0.5,0.5) 2.718348 2.718282 6.575685e-05 240 (0.5,0.6) 3.004239 3.004166 7.289348e-05 280 (0.5,0.7) 3.320198 3.320117 8.064224e-05 320 (0.5,0.8) 3.669386 3.669297 8.915426e-05 360 (0.5,0.9) 4.055299 4.055200 9.854220e-05 400 (0.5,1.0) 4.481798 4.481689 1.089103e-04③取不同步长时数值解的最大误差④误差曲线图⑤误差曲面图图中：上曲面h=1/10 τ=1/100下曲面h=1/20 τ=1/400（3）MATLAB源代码h=input('输入空间步长h=');tau=input('输入时间步长tau=');r=tau/h^2;m=1/h;n=1/tau;i=1:m-1;k=1:n-1;A=diag(-r*ones(m-2,1),-1)+diag((1+2*r)*ones(m-1,1))+diag(-r*ones(m-2,1),1);B=diag(r*ones(m-2,1),-1)+diag((1-2*r)*ones(m-1,1))+diag(r*ones(m-2,1),1);xi=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;ui0=exp(xi)';u00=1;%u01=1+tau;um0=exp(1);u0k=exp(tk)';umk=exp(1+tk)';fik=zeros(m-1,1);uik=zeros(m-1,n);%uik2=ones(m-1,n-1);uik(:,1)=(1+tau)*exp(xi);%first k=1%C=B*ui0+2*tau*fik+[r*(u0k(2)+r*u00);zeros(m-3,1);r*(umk(2)+um0)];%uik(:,2)=tridiag(A,C,m-1);uik(:,2)=inv(A)*C;%%%%%%%%%%%%%%%%for k=2:n-1C=B*uik(:,k-1)+2*tau*fik+[r*(u0k(k+1)+u0k(k-1));zeros(m-3,1);r*(umk(k+1)+umk(k-1))];uik(:,k+1)=inv(A)*C;endU2=exp(xi+1)';true=zeros(m-1,n);for i=1:m-1true(i,:)=exp(xi(i)+tk);endfprintf('h=%.3f tau=%.3f\n',h,tau);fprintf('k (x,t) 数值解精确解误差\n');for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));elsefprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));endenda=max(abs(uik-true));E=max(a);fprintf('Emax=%e',E);% e=max(true-uik); % h=max(e)%plot(e)%plot(abs(uik(:,n)-true(:,n)))第七题补充题2 非线性方程（教材P129/算例3.4）（1）差分格式ðu ðt −u4ð2uðt2=f(x,t)u i k+1+u i k−12τ−(u i k)4δx2(u i k−1+u i k+12)=f(x i,t k)（2）取不同步长时数值解的最大误差（3）误差曲线图（4）误差曲面图图中：上曲面h=τ=1/10中曲面h=τ=1/20下曲面h=τ=1/40（5）MATLAB源代码h=input('输入空间步长h=');tau=input('输入时间步长tau=');r=tau/h^2;m=1/h;n=1/tau;i=1:m-1;k=1:n-1;xi=0+1/m:1/m:1-1/m;tk=1/n:1/n:1;ui0=sin(pi/4+xi*pi/2)';u00=sqrt(2)/2;um0=sqrt(2)/2;u0k=sqrt(2)/2*exp(-tk)';umk=sqrt(2)/2*exp(-tk)';fik=zeros(m-1,n);for j=1:nfik(:,j)=(exp(-tk(j))*sin(pi/2+xi*pi/2).*(pi^2/4*exp(-4*tk(j)).*sin(pi/4+xi*pi/2).^4-1))'; enduik=zeros(m-1,n);%uik2=ones(m-1,n-1);uik(:,1)=(sin(pi/4+xi*pi/2)+tau*(sin(pi/4+xi*pi/2).^4).*(-pi^2/4*sin(pi/4+xi*pi/2))+tau*sin(pi/4+xi*pi/2).*(pi^2/4*sin(pi/4+xi*pi/2).^4-1))';%first k=1%A=zeros(m-1,m-1);B=zeros(m-1,m-1);A(1,1)=1+2*r*(uik(1,1)^4);A(1,2)=-r*(uik(1,1)^4);B(1,1)=1-2*r*(uik(1,1)^4);B(1,2)=r*(uik(1,1)^4);A(m-1,m-2)=-r*(uik(m-1,1)^4);A(m-1,m-1)=1+2*r*(uik(m-1,1)^4);B(m-1,m-2)=r*(uik(m-1,1)^4);B(m-1,m-1)=1-2*r*(uik(m-1,1)^4);for i=2:m-2A(i,i-1)=-r*(uik(i,1)^4);A(i,i)=1+2*r*(uik(i,1)^4);A(i,i+1)=-r*(uik(i,1)^4);B(i,i-1)=r*(uik(i,1)^4);B(i,i)=1-2*r*(uik(i,1)^4);B(i,i+1)=r*(uik(i,1)^4);endC=B*ui0+2*tau*fik(:,1)+[r*(uik(1,1)^4*(u0k(2)+u00));zeros(m-3,1);r*(uik(m-1,1))^4*(umk(2)+um0)];%uik(:,2)=tridiag(A,C,m-1);uik(:,2)=inv(A)*C;%%%%%%%%%%%%%%%%for k=2:n-1A=zeros(m-1,m-1);B=zeros(m-1,m-1);A(1,1)=1+2*r*(uik(1,k)^4);A(1,2)=-r*(uik(1,k)^4);B(1,1)=1-2*r*(uik(1,k)^4);B(1,2)=r*(uik(1,k)^4);A(m-1,m-2)=-r*(uik(m-1,k)^4);A(m-1,m-1)=1+2*r*(uik(m-1,k)^4);B(m-1,m-2)=r*(uik(m-1,k)^4);B(m-1,m-1)=1-2*r*(uik(m-1,k)^4);for i=2:m-2A(i,i-1)=-r*(uik(i,k)^4);A(i,i)=1+2*r*(uik(i,k)^4);A(i,i+1)=-r*(uik(i,k)^4);B(i,i-1)=r*(uik(i,k)^4);B(i,i)=1-2*r*(uik(i,k)^4);B(i,i+1)=r*(uik(i,k)^4);endC=B*uik(:,k-1)+2*tau*fik(:,k)+[r*(uik(1,k)^4)*(u0k(k+1)+u0k(k-1));zeros(m-3,1);r*(uik(m-1,k)^4)*(umk(k+1)+umk(k-1))];uik(:,k+1)=inv(A)*C;endtrue=zeros(m-1,n);for i=1:m-1true(i,:)=exp(-tk).*(sin(pi/4+xi(i)*pi/2));end% e=max(true-uik);% h=max(e)% plot(e)% plot(uik(:,90));% hold on;plot(true(:,end)-uik(:,end))%surf(abs(uik-true))fprintf('h=%.3f tau=%.3f\n',h,tau);fprintf('k (x,t) 数值解精确解误差\n'); for k=1:10p=int16(k*0.1/tau);if p<100fprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));elsefprintf('%d(0.5,%.1f) %f %f %e\n',p,0.1*k,uik(m/2,p),true(m/2,p),abs(uik(m/2,p)-true(m/2,p)));endenda=max(abs(uik-true));E=max(a);fprintf('Emax=%e',E);。
