孙志忠北京理工大学偏微分方程数值解上机作业
姓名_学号_中国海洋大学偏微分方程的数值解法第七次作业
偏微分方程的数值解法上机习题七题目：求解边值问题：()()()()()()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 超松弛法最少。
北理研究生数值分析---第九章课件
§1 Euler方法
1.Euler方法 以差分方程初值问题
y n 1 y n hf ( x n , y n ) y 0 y (a ) n 0 ,1, , N 1
的解作为微分方程初值问题的数值解，即
y(xn ) yn
称为Euler方法。
1.用Euler方法求初值问题
( k 0 ,1, 2 , )
y n 1 y n 1
(k )
( k 1)
h 2
f ( x n 1 , y n 1 ) f ( x n 1 , y n 1 )
(k )
( k 1 )
hL 2 hL 2
y n 1 )
k
( k 1)
y n 1
(k ) (0)
k
hL 2
( )
k p 1
hL 2
)
k
] y n 1 y n 1
(1 ) (0)
1
hL 2
y n 1 y n 1 0
(1 ) (0)
2.2 改进Euler法
y n 1 y n hf ( x n , y n ) h y n 1 y n f ( x n , y n ) f ( x n 1 , y n 1 ) 2 预测 校正
约定：不加特别说明，必有 （1） f ( x , y ) 连续，且关于 y 满足 Lipschitz 条件，由此保证初始问题的解存在唯一。 （2）步长 h n x n 1 x n ( n 0 ,1, , N 1 ) 为常量 h 。
微分方程离散化的方法
(1) 用差商近似导数
y ( x n ) y ( x n 1 ) y ( x n ) h
(完整word版)偏微分方程数值解习题解答案
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）的解。
2016-偏微分方程数值解法-课程大纲-谢树森
中国海洋大学本科生课程大纲课程属性：公共基础/通识教育/学科基础/专业知识/工作技能，课程性质：必修、选修一、课程介绍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，这样可使边界处的结果也为中心差分，保证了格式的一致精度。
偏微分方程数值解上机实验报告(matlab做的)
偏微分方程数值解法上机报告（一）一、实验题目：用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. 编程并作图，将计算解与精确解比较。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
偏微分方程数值解大作业目录第一题 (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);。