偏微分方程数值及matlab实验报告.docx
偏微分方程数值解第三次上机实验报告
偏微分方程数值解法第三次上机实验报告一、实验题目：用线性元求解下列问题的数值解：2,1,1(x,1)u(x,1)0,1x 1(1,y)1,u (1,y)0,11xx u x y u u y =--<<⎧⎪-==-<<⎨⎪-==-<<⎩ （精确到小数点后四位）二、实验过程：利用PDEToolbox 工具箱求解该偏微分方程。
分析：方程是Possion 方程形式c u au f -+=，其中c=-1，a=0，f=-2； 第一个边界条件是Dirichlet 条件，第二个边界条件是Neumann 条件。
1.在MA TLAB 命令窗口键入pdetool 并运行，打开PDEToolbox 界面；2.在Options 菜单下选择Grid 命令，显示网格，能更容易确定所绘图形的大小；3.绘出区域，选择Boundary 的Boundary Mode ，双击每个边界，设置边界相应的参数值；4.选择PDE 菜单中PDE Mode 命令，进入PDE 模式。
单击PDE 菜单中PDE Specification ….选项，设置方程类型及参数；5.选择Mesh 菜单中Initialize Mesh 命令，进行网格剖分，再选择Refine Mesh 命令，进行网格加密，如下图：三、实验结果：选择Solve 菜单中solve PDE 命令，解偏微分方程，其图形解如图：图1 图形解图2 三维图形解图3 解的等值线图和矢量场图选择Mesh菜单中的Export Mesh，得到结点xy坐标；选择Solve菜单中的Export Solution…，得到每个节点处的值，输出u，即解的数值。
四、实验总结：通过本次试验，掌握了利用Matlab中的PDE求解工具得到PDE的解的方法，并对偏微分方程的形式有了更多的掌握。
用matlab对微分方程求解实验报告.
o 《高等数学》上机作业(三一、上机目的1、学会用 M a t l a b 求简单微分方程的解析解。
2、学会用 M a t l a b 求微分方程的数值解。
二、上机内容1、求简单微分方程的解析解.2、求微分方程的数值解.3、数学建模实例.4、上机作业. 三、上机作业1. 求微分方程:在初值条件下的特解,并画出解函数的图形. 命令>> y =d s o l v e ('x *D y +y -e x p (x =0','y (1=2*e x p (1','x ' 运行结果:y = 1/x *e x p (x +1/x *e x p (1'xxy y e +-=12(y e =函数图象:2. 求微分方程的特解.22450(00,'(110d y dyy dx dx y y ?+-=???==?命令>> y=dsolve('D2y+4*Dy-5*y=0','y(0=0,Dy(1=10','x' 运行结果:y=10/(exp(1+5*exp(-5*exp(x-10/(exp(1+5*exp(-5*exp(-5*x3. 鱼雷追击问题一敌舰在某海域内沿着正北方向航行时,我方战舰恰好位于敌舰的正西方向 1 公里处.我舰向敌舰发射制导鱼雷,敌舰速度为0.42 公里/分,鱼雷速度为敌舰速度的2倍。
试问敌舰航行多远时将被击中?M文件x0=0; xf=0.9999999999999; [x,y]=ode15s('eq1',[x0 xf],[0 0]; plot(x,y(:,1,'b.'hold on;y=0:0.1:1;plot(1,y, '*'运行结果图像:结论:大概在y=0.67处击中敌方舰艇!(选做一个慢跑者在平面上沿椭圆以恒定的速率v=1跑步,设椭圆方程为:x=10+20cost, y=20+5sint. 突然有一只狗攻击他. 这只狗从原点出发,以恒定速率w跑向慢跑者,狗的运动方向始终指向慢跑者.分别求出w=20,w=5时狗的运动轨迹.W=20M文件代码function dy=eq3(t,ydy=zeros(2,1;dy(1=20*(10+20*cos(t-y(1/sqrt((10+20*cos(t-y(1^2+(2 0+15*sin(t-y(2^2;dy(2=20*(20+15*sin(t-y(2/sqrt((10+20*cos(t-y(1^2+(2 0+15*sin(t-y(2^2;运行命令t0=0;tf=10;[t,y]=ode45('eq3',[t0 tf],[0 0];T=0:0.1:2*pi;X=10+20*cos(T;Y=20+15*sin(T;plot(X,Y,'-'hold onplot(y(:,1,y(:,2,'r*'运行结果:利用二分法更改tf tf=5时tf=2.5时tf=3.15时:所以在t=3.15时刻恰好追上!W=5M文件代码function dy=eq4(t,ydy=zeros(2,1;dy(1=5*(10+20*cos(t-y(1/sqrt((10+20*cos(t-y(1^2+(20 +15*sin(t-y(2^2; dy(2=5*(20+15*sin(t-y(2/sqrt((10+20*cos(t-y(1^2+(20 +15*sin(t-y(2^2; 命令:t0=0;tf=10;[t,y]=ode45('eq4',[t0 tf],[0 0]; T=0:0.1:2*pi;X=10+20*cos(T; Y=20+15*sin(T; plot(X,Y,'-'hold onplot(y(:,1,y(:,2,'*' 运行结果更改tf=20运行结果Tf=40 11所以永远追不上！四、上机心得体会高等数学是工科学生的主干科目，它应用于生产生活的方方面面，通过建模，计算可以求出实际问题的最优化问题！因此我们需要掌握建模和利用专业软件处理实际问题的能力！ 12。
偏微分方程数值及matlab实验报告
偏微分方程数值实验报告八实验题目：利用有限差分法求解.0)1(,0)1(),()()(==-=+''-u u x f x u x u 真解为)1()(22x e x u x -=-实现算法：对于两点边值问题,)(,)(,,d 22βα==∈=-b u a u l x f dxu(1)其中),(b a l =f b a ),(<为],[b a l =上的连续函数，βα，为给定常数.其相应的有限差分法的算法如下：1.对求解区域做网格剖分，得到计算网格.在这里我们对区间l 均匀剖分n 段，每个剖分单元的剖分步长记为nab h -=.2.对微分方程中的各阶导数进行差分离散，得到差分方程.运用的离散方法有：方法一:用待定系数和泰勒展开进行离散)()()()(d )(d 111122++--++≈i i i i i i i i x u x u x u x x u ααα方法二：利用差商逼近导数21122)()(2)()(d )(d h x u x u x u x x u i i i i i -++-≈(2)将(2)带入(1)可以得到)()()()(2)(211u R x f h x u x u x u i i i i i +=+---+,其中)(u R i 为无穷小量，这时我们丢弃)(u R i ，则有在i x 处满足的计算公式：1,...,1)()()(2)(211-==+---+n i x f hx u x u x u i i i i ，(3)3.根据边界条件，进行边界处理.由(1)可得βα==n u u ,0(4)称(3)(4)为逼近(1)的差分方程，并称相应的数值解向量1-n U 为差分解，i u 为)(i x u 的近似值.4.最后求解线性代数方程组，得到数值解向量1-n U .实验题目：用Lax-Wendroff 格式求解方程：.4sin 1),1(],1,0[,2sin 1)0,(,0),1,0(,02t t u x x x u t x u u x t ππ+=∈+=>∈=- （1） （精确解).2(2sin 1t x u ++=π） 数值边值条件分别为： （a ）).(20101n 0nn nu u hu u -+=+τ （b ）.1n 0nu u =（c ）.02-12111n 0=++++n n u u u请将计算结果与精确解进行比较。
二阶椭圆偏微分方程实例求解(附matlab代码)
《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名：学号：班级：2013年11月19日二阶椭圆偏微分方程第一边值问题摘要对于解二阶椭圆偏微分方程第一边值问题.课本上已经给出了相应的差分方程。
而留给我的难题就是把差分方程组表示成系数矩阵的形式.以及对系数进行赋值。
解决完这个问题之后.我在利用matlab 解线性方程组时.又出现“out of memory ”的问题。
因为99*99阶的矩阵太大.超出了分配给matlab 的使用内存。
退而求其次.当n=10.h=1/10或n=70.h=1/70时.我都得出了很好的计算结果。
然而在解线性方程组时.无论是LU 分解法或高斯消去法.还是gauseidel 迭代法.都能达到很高的精度。
关键字：二阶椭圆偏微分方程 差分方程 out of memory LU 分解 高斯消去法 gauseidel 迭代法一、题目重述解微分方程：()()2222((,))((,))()(,)()(,)(,)1y x x x y y x y yxxyxye u x y e u x y x y u x y x y u x y u x y y e x e e y x e--+++-+=-++++已知边界:(0,)1,(1,),(,0)1,(,1)y x u y u y e u x u x e ====求数值解, 把区域[0,1][0,1]G =?分成121/100,1/100h h ==.n =100 注：老师你给的题F 好像写错了.应该把22x y y e x e +改成22y x y e x e +。
二、问题分析与模型建立2.1微分方程上的符号说明()()22221y x xy xy y e x e e y x e -++++2.2课本上差分方程的缺陷课本上的差分方程为：举一个例子：当i=2,j=3时.；当i=3,j=3时.。
但是.显然这两个不是同一个数.其大小也不相等。
matlab求微分方程的解 实验报告四
对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行高中资料试卷调整试验；通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术关，通系电1，力过根保管据护线生高0不产中仅工资2艺料22高试2可中卷以资配解料置决试技吊卷术顶要是层求指配，机置对组不电在规气进范设行高备继中进电资行保料空护试载高卷与中问带资题负料2荷试2，下卷而高总且中体可资配保料置障试时2卷，32调需3各控要类试在管验最路；大习对限题设度到备内位进来。行确在调保管整机路使组敷其高设在中过正资程常料1工试中况卷，下安要与全加过，强度并看工且25作尽52下可22都能护可地1关以缩于正小管常故路工障高作高中；中资对资料于料试继试卷电卷连保破接护坏管进范口行围处整，理核或高对者中定对资值某料，些试审异卷核常弯与高扁校中度对资固图料定纸试盒，卷位编工置写况.复进保杂行护设自层备动防与处腐装理跨置，接高尤地中其线资要弯料避曲试免半卷错径调误标试高方中等案资，，料要编试求5写、卷技重电保术要气护交设设装底备备置。4高调、动管中试电作线资高气，敷料中课并设3试资件且、技卷料中拒管术试试调绝路中验卷试动敷包方技作设含案术，技线以来术槽及避、系免管统不架启必等动要多方高项案中方；资式对料，整试为套卷解启突决动然高过停中程机语中。文高因电中此气资，课料电件试力中卷高管电中壁气资薄设料、备试接进卷口行保不调护严试装等工置问作调题并试，且技合进术理行，利过要用关求管运电线行力敷高保设中护技资装术料置。试做线卷到缆技准敷术确设指灵原导活则。。：对对在于于分调差线试动盒过保处程护，中装当高置不中高同资中电料资压试料回卷试路技卷交术调叉问试时题技，，术应作是采为指用调发金试电属人机隔员一板，变进需压行要器隔在组开事在处前发理掌生；握内同图部一纸故线资障槽料时内、，设需强备要电制进回造行路厂外须家部同出电时具源切高高断中中习资资题料料电试试源卷卷，试切线验除缆报从敷告而设与采完相用毕关高，技中要术资进资料行料试检，卷查并主和且要检了保测解护处现装理场置。设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况，然后根据规范与规程规定，制定设备调试高中资料试卷方案。
matlab求微分方程的解-实验报告四
matlab求微分方程的解-实验报告四《matlab与数学实验》实验报告实验序号：实验四日期： 2015年 5 月 25 日班级132132002姓名彭婉婷学号 1321320057 实验名称求微分方程的解问题背景描述实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．实验目的本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍 Euler 折线法．实验原理与数学模型MATLAB7.11.0主要内容（要点）1. 求微分方程0sin2')1(2=-+-xxyyx的通解．2. 求微分方程xeyyy x sin5'2''=+-的通解．3. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=-+=++yxdtdyyxdtdx在初始条件0|,1|====ttyx下的特解，并画出解函数()y f x=的图形．4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t∈．利用画图来比较两种求解器之间的差异．5. 用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=-=1)0(,12'32yyxyy的数值解(步长h取0.1)，求解范围为区间[0,2]．选做：6. 用四阶 Runge-Kutta 法求解微分方程初值问题⎩⎨⎧=-=1)0(,cos'yxeyy x的数值解(步长h取0.1)，求解范围为区间[0,3]．迭代法实验过程记录（含基本步骤、主要程序清单及异常情况记录等）1.求微分方程0sin2')1(2=-+-xxyyx的通解．程序:clearsyms x yy=dsolve('(x^2-1)*Dy+2*x*y=sin(x)','x') 答案:y =-(C2 + cos(x))/(x^2 - 1)2.求微分方程xeyyy x sin5'2''=+-的通解．程序:clearsyms x yy=dsolve('D2y-2*Dy+5*y=exp(x)*sin(x)','x ')simplyify(x/y)weijiao答案:y =(exp(x)*sin(x))/6 - (sin(3*x)*exp(x))/8 + (sin(5*x)*exp(x))/24 + C4*cos(2*x)*exp(x) + sin(2*x)*exp(x)*(cos(x)/4 - cos(3*x)/12 + 1/6) + C5*sin(2*x)*exp(x)3. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dt dy y x dt dx在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =的图形．程序:clearsyms x y t[x,y]=dsolve('Dx+x+y=0','Dy+x-y=0','x(0)=1','y(0)=0','t')ezplot(x,y,[0, 1])(t 的取值，t 是与x,y 相关的，如果不给范围，则会默认为一个较大的区间) simplify(x)simplify(y)答案:x =exp(2^(1/2)*t)/2 + 1/(2*exp(2^(1/2)*t)) - (2^(1/2)*exp(2^(1/2)*t))/4 +2^(1/2)/(4*exp(2^(1/2)*t))y =2^(1/2)/(4*exp(2^(1/2)*t)) -(2^(1/2)*exp(2^(1/2)*t))/4图形:4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为t ．利用画图来比较两种求解器之间的差异．[0,2]先编写函数文件verderpol.m：clearfunction xprime=verderpol(t,x)xprime=[-x(1)-x(2); x(2)-x(1)];再编写命令文件cleary0=[1;0];[t,x] = ode45('verderpol',[0,2],y0);x1=x(:,1);x2=x(:,2);plot(x1,x2,'b-')hold ony0=[1;0];[t,x]=ode23('verderpol',[0,2],y0);x1=x(:,1);x2=x(:,2);plot(x(:,1),x(:,2),'r-');图形:两种求解器之间的差异:ode45大部分场合的首选算法ode23使用于精度较低的情形但在此题中并没有体现出差异。
偏微分方程(PDEs)的MATLAB数值解法
偏微分方程的MATLAB求解精讲©MA TLAB求解微分/偏微分方程，一直是一个头大的问题，两个字，“难过”，由于MA TLAB对LaTeX的支持有限，所有方程必须化成MA TLAB可接受的标准形式，不支持像其他三个数学软件那样直接傻瓜式输入，这个真把人给累坏了！不抱怨了，还是言归正传，回归我们今天的主体吧！MA TLAB提供了两种方法解决PDE问题，一是pdepe()函数，它可以求解一般的PDEs，据用较大的通用性，但只支持命令行形式调用。
二是PDE工具箱，可以求解特殊PDE问题，PDEtool有较大的局限性，比如只能求解二阶PDE问题，并且不能解决偏微分方程组，但是它提供了GUI界面，从繁杂的编程中解脱出来了，同时还可以通过File->Save As直接生成M代码一、一般偏微分方程组(PDEs)的MA TLAB求解 (3)1、pdepe函数说明 (3)2、实例讲解 (4)二、PDEtool求解特殊PDE问题 (6)1、典型偏微分方程的描述 (6)（1）椭圆型 (6)（2）抛物线型 (6)（3）双曲线型 (6)（4）特征值型 (7)2、偏微分方程边界条件的描述 (8)（1）Dirichlet条件 (8)（2）Neumann条件 (8)3、求解实例 (9)一、一般偏微分方程组(PDEs)的MATLAB 求解1、pdepe 函数说明MA TLAB 语言提供了pdepe()函数，可以直接求解一般偏微分方程(组)，它的调用格式为sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)【输入参数】@pdefun ：是PDE 的问题描述函数，它必须换成下面的标准形式(,,)[(,,,)](,,,)()m m u u u uc x t x x f x t u s x t u x t x x x−∂∂∂∂∂=+∂∂∂∂∂式1 这样，PDE 就可以编写下面的入口函数 [c,f,s]=pdefun(x,t,u,du)m,x,t 就是对应于(式1)中相关参数，du 是u 的一阶导数，由给定的输入变量即可表示出出c,f,s 这三个函数@pdebc ：是PDE 的边界条件描述函数，必须先化为下面的形式(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂+=∂ 于是边值条件可以编写下面函数描述为 [pa,qa,pb,qb]=pdebc(x,t,u,du)其中a 表示下边界，b 表示下边界@pdeic ：是PDE 的初值条件，必须化为下面的形式00(,)u x t u =我们使用下面的简单的函数来描述为 u0=pdeic(x)m,x,t ：就是对应于(式1)中相关参数【输出参数】sol ：是一个三维数组，sol(:,:,i)表示u i 的解，换句话说u k 对应x(i)和t(j)时的解为sol(i,j,k)通过sol ，我们可以使用pdeval()直接计算某个点的函数值2、实例讲解试求解下面的偏微分2111222221220.024()0.17()u u F u u t xu u F u u tx ∂∂=−− ∂∂ ∂∂ =−− ∂∂ 其中， 5.7311.46()x x F x e e −=−，且满足初始条件12(,0)1,(,0)0u x u x ==及边界条件1221(0,)0,(0,)0,(1,)1,(1,)0u ut u t u t t x x∂∂====∂∂【解】（1）对照给出的偏微分方程，根据标注形式，则原方程可以改写为111222120.024()1.*1()0.17u u F u u x u u F u u t t x ∂−−∂∂∂=+ ∂−∂∂∂可见m=0，且1122120.024()1,,1()0.17u F u u x c f s u F u u x ∂−− ∂===∂−∂%% 目标PDE 函数function [c,f,s]=pdefun (x,t,u,du) c=[1;1];f=[0.024*du(1);0.17*du(2)]; temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp));（2）边界条件改写为12011010.*.*00000u f f u −+=+=下边界上边界%% 边界条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t) %a 表示下边界，b 表示上边界 pa=[0;ua(2)];qa=[1;0]; pb=[ub(1)-1;0]; qb=[0;1];（3）初值条件改写为1210u u =%% 初值条件函数function u0=pdeic(x) u0=[1;0];（4）最后编写主调函数 clcx=0:0.05:1; t=0:0.05:2; m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);figure('numbertitle','off','name','PDE Demo ——by Matlabsky') subplot(211)surf(x,t,sol(:,:,1)) title('The Solution of u_1') xlabel('X') ylabel('T') zlabel('U') subplot(212)surf(x,t,sol(:,:,2)) title('The Solution of u_2') xlabel('X') ylabel('T') zlabel('U')二、PDEtool 求解特殊PDE 问题MATLAB 的偏微分工具箱(PDE toolbox)可以比较规范的求解各种常见的二阶偏微分方程，但是惋惜的是只能求解特殊二阶的PDE 问题，并且不支持偏微分方程组！PDE toolbox 支持命令行形式求解PDE 问题，但是要记住那些命令以及调用形式真的很累人，还好MATLAB 提供了GUI 可视交互界面pdetool ，在pdetool 中可以很方便的求解一个PDE 问题，并且可以帮我们直接生成M 代码(File->Save As)。
(完整word版)Matlab数学实验报告
Matlab 数学实验报告一、实验目的通过以下四组实验,熟悉MATLAB的编程技巧,学会运用MATLAB的一些主要功能、命令，通过建立数学模型解决理论或实际问题。
了解诸如分岔、混沌等概念、学会建立Malthu模型和Logistic 模型、懂得最小二乘法、线性规划等基本思想。
二、实验内容2.1实验题目一2.1.1实验问题Feigenbaum曾对超越函数y=λsin(πｘ)（λ为非负实数）进行了分岔与混沌的研究，试进行迭代格式x k+1=λsin(πx k),做出相应的Feigenbaum图2.1.2程序设计clear;clf;axis([0,4,0,4]);hold onfor r=0:0.3:3.9x=[0.1];for i=2:150x(i)=r*sin(3.14*x(i-1));endpause(0.5)for i=101:150plot(r,x(i),'k.');endtext(r-0.1,max(x(101:150))+0.05,['\it{r}=',num2str(r)]) end加密迭代后clear;clf;axis([0,4,0,4]);hold onfor r=0:0.005:3.9x=[0.1];for i=2:150x(i)=r*sin(3.14*x(i-1));endpause(0.1)for i=101:150plot(r,x(i),'k.');endend运行后得到Feigenbaum图2.2实验题目二2.2.1实验问题某农夫有一个半径10米的圆形牛栏，长满了草。
他要将一头牛拴在牛栏边界的桩栏上，但只让牛吃到一半草，问拴牛鼻子的绳子应为多长？2.2.2问题分析如图所示，E为圆ABD的圆心，AB为拴牛的绳子，圆ABD为草场，区域ABCD为牛能到达的区域。
问题要求区域ABCD等于圆ABC的一半，可以设BC等于x,只要求出∠a和∠b就能求出所求面积。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
偏微分方程数值实验报告八实验题目：利用有限差分法求解u ( x) u(x) f (x),u( 1) 0, u(1)0.真解为u( x)e x2(1 x 2)实现算法：对于两点边值问题d 2u f , xl ,dx 2(1)u(a),u(b),其中 l ( a, b) (ab), f 为 l[ a,b] 上的连续函数， ，为给定常数 .其相应的有限差分法的算法如下：1.对求解区域做网格剖分，得到计算网格.在这里我们对区间 l 均匀剖分 n 段，每个剖分单元ba的剖分步长记为 h .n2.对微分方程中的各阶导数进行差分离散，得到差分方程.运用的离散方法有：方法一 :用待定系数和泰勒展开进行离散d 2u( x i ) i 1u( x i 1)i u( x i )i 1u( x i 1)d( x i )2方法二：利用差商逼近导数d 2u( x i ) u( x i 1 ) 2u( x i ) u( x i 1 ) d( x i )2 h 2将(2) 带入 (1)可以得到u(x i 1) 2u(x i ) u(x i 1 )) R i (u) ,h 2f ( x i其中 R i (u) 为无穷小量，这时我们丢弃 R i (u) ，则有在 x i 处满足的计算公式：u(x i 1) 2u( x i ) u( x i 1 )1,..., n 1h 2 f ( x i )， i3.根据边界条件，进行边界处理.由 (1)可得u 0, u n(2)(3)(4)称(3)(4)为逼近 (1) 的差分方程，并称相应的数值解向量U n 1 为差分解， u i 为 u( x i ) 的近似值 .4.最后求解线性代数方程组，得到数值解向量U n 1 .实验题目： 用 Lax-Wendroff格式求解方程：u t 2u x 0, x (0,1), t 0,u( x,0) 1 sin 2 x, x [ 0,1],u(1,t ) 1 sin 4 t.（精确解 u 1 sin 2 ( x2t). ）数值边值条件分别为：（a ） u 0n 1u 0n2(u 1n u 0n ).h（b ） u 0n u 1n .（c ） u 0n 1 - 2u 1n 1 u 2n 1 0.请将计算结果与精确解进行比较。
实现算法： 1. 网格剖分：对求解区域 G (0,1) (0,T ] 作均匀网格剖分 .节点：x jjh , j 0,1,..., Nt k jh ,k0,1,..., M其中空间和时间步长：h1 , T .NM2. 算法实现将 u(x i , t k 1 ) 在节点 ( x i ,t k ) 处作泰勒级数展开u(x i , t k 1) u(x i , t k )[u k 22u kO ( 3)]i[t 2 ]it2!考虑在节点 ( x j , t k ) 处（ 1）的微分方程，有：u a u 0.t x2u( a u a ( u22ut 2) ) ax 2.txx t将上述两式代入（ 2）式，得u( x i ,t k 1 ) u( x i ,t k ) a [ u] i k a 2 2 [ 2 u 2 ]i k O ( 3 )x 2 x对 x 的一阶、二阶导数用中心差商代替（ 1）（ 2）u k12[x ]i2h[( u( x i1, t k )u( x i 1, t k ))]O(h)[2u k12u( x i ,t k )u( x i 1, t k ))]2) x2 ]i h2 [( u(x i 1 ,t k )O (h代入整理后得到u( x i ,t k 1 ) u( x i ,t k )a[ u( x i 1 ,t k ) u( x i 1 ,t k )] 2ha 22O( h2 ) 2 h2 )O( 3 ) 2[ u( x i 1 ,t k )2u( x i , t k )u( x i1 , t k )]O (2h略去误差项，以u i k代替 u( x i ,t k ) ，得到如下差分格式a(u i k122u i k 1u i k u i k1 )a2 (u i k12u i k u i k1)（ 3）2h2h（ 3）式就是 Lax-Wendroff 格式，其截断误差为O(2h2 ) ，节点如图| a |，就得到（1）式的 Lax-Wendroff格式的公式令 r hu i k 1u i k r(u i k1u i k1 )r 2(u i k12u i k u i k1)（ 4）22（ 4）式是二阶精度的差分格式.程序代码：function [X,T,U] = advection_fd1d (NS ,NT ,pde,bd)%WAVE_EQUATION_FD1D利用有限差分方法计算一维双曲线方程%输入参数：%NS整型，空间剖分段数%NT整型，时间剖分段数%pde结构体，带求解的微分方程模型的已知数据，%如边界、初始、系数和右端项等条件. %bd数值边值条件%输出参数 :%X长度NS+1的列向量，空间网格剖分%T长度NT+1的行向量，时间网格剖分%U (NS+1)*(NT+1)矩阵，U(:,i)表示第i个时间层网格剖分上的数值解[X,h] = (NS);[T,tau] = (NT);N = length(X); M = length(T);U = zeros(N,M);%初值条件U(:,1) = (X);a = ;r = a*tau/h;%边值条件if a >= 0% 左边值条件U(1,:) = (T)elseU(end,:) = (T)%右边值条件endfor i = 2:MU(2:end -1,i) =U(2:end-1,i-1)-r*(U(3:end,i-1)-U(1:end-2,i-1))/2+...r^2*(U(3:end,i-1)-2*U(2:end-1,i-1)+U(1:end-2,i-1))/2;switch (bd)case { 'a0' }a0();case { 'b' }b();case { 'c' }c();otherwisedisp(['Sorry, I do not know your ', bd]);endendfunction a0()U(1,i)=U(1,i-1)-r*(U(2,i-1)-U(1,i-1));endfunction b()U(1,i)=U(2,i-1);endfunction c()U(1,i)=2*U(2,i)-U(3,i);endendfunction pde = model_data () %MODEL_DATA数据模型TI = 0;TF = 1;SI = 0;SF = 1;pde = struct('u_exact' ,@u_exact, 'u_initial' 'u_left',@u_left, 'u_right',@u_right,,@u_initial,'time_grid', ......@time_grid,'space_grid',@space_grid, 'advection_fd1d_error',@advection_fd1d_error,'a ' ,-2);function [T,tau] = time_grid(NT)T = linspace(TI,TF,NT+1);tau = (TF-TI)/NT;endfunction [X,h] = space_grid(NS)X = linspace(SI,SF,NS+1)'h = (SF-SI)/NS;endfunction U = u_exact(X,T)[x,t]=meshgrid(X,T);U = 1+sin(2*pi*(x+2*t));endfunction u = u_initial (x)u = 1+sin(2*pi*x);endfunction u = u_right(t)u=1+sin(4*pi*t);endendfunction showsolution (X,T,U)%% SHOWSOLUTION以二元函数方式显示数值解% 输入参数% X 长度为 NS +1的列向量，空间网格剖分N% T 长度为 NT +1的行向量，时间网格剖分M%U N*M 矩阵， U(:,i) 表示第 i 个时间层网格部分上的数值解[x,t] = meshgrid (X,T);mesh (x,t,U');xlabel ( 'X' );ylabel ( 'T' );zlabel ( 'U(X,T)' );endfunction showvarysolution (X,T,U,UE)%% SHOWVARYSOLUTION显示数值解随着时间的变化%输入参数% X长度为NS +1的列向量，空间网格剖分N% T长度为NT +1的行向量，时间网格剖分M%U N*M矩阵，U(:,i)表示第i个时间层网格部分上的数值解M = size (U ,2) ;figurexlabel ( 'X' );ylabel ( 'U' );s = [X(1) ,X(end),min(min(U)),max(max(U))];axis (s);for i = 1:Mplot (X,U(:,i));axis (s);pause ;title (['T=' ,num2str(T(i)),' 时刻的温度分布' ]) End%一维双曲线有限差分方法主测试脚本pde=model_data()[X,T,U]=advection_fd1d(100,200,pde,'a' );UE=(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解[X,T,U]=advection_fd1d(100,200,pde,'b' );UE=(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解'c' );[X,T,U]=advection_fd1d(100,200,pde,UE=(X,T);showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解showsolution(X,T,U);%以二元函数方式显示数值解。