偏微分上机报告一

偏微分方程数值解法第三次上机实验报告一、实验题目：用线性元求解下列问题的数值解：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的解的方法，并对偏微分方程的形式有了更多的掌握。

偏微分方程数值解大作业目录第一题 (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)可作出误差曲线注意精度不要选取的太低，否则随着步长减小误差反而增大。

偏微分方程数值解上机实习数值求解二维扩散方程的初边值问题2222(,,0)sin sin 2(0,,)(1,,)1(,0,)(,1,)1u u ut x yu x y x y u y t u y t u x t u x t ππ⎧∂∂∂=+⎪∂∂∂⎪⎪=⎨⎪==⎪⎪==⎩(0,1,0)(0,1)(01,0)(01,0)x y t x y y t x t <<><<<<≥<<≥ 古典显式格式：1,,1,,1,,1,,12222n nn n nn n nj l j lj l j l j lj l j l j l u u u u u u u u hhτ++-+---+-+=+将原格式化为：12,1,1,,1,1,(14),/n nnn nnj l j l j l j l j l j l u u u u u u h λλλλλλτ++-+-=++++-=其中 附源程序：%-------------------------------------------运用古典显式差分格式求解二维扩散方程的初边值问题; function gdxs(ti,h,t)%-------------------------------------------ti:时间步长; %-------------------------------------------h:空间步长; k=t/ti;m=1/h+1;r=ti/h^2; %------------------------------ r 为网格比; w=ones(m,m); u=ones(m,m); for i=2:m-1for j=2:m-1u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1)); end end ticfor l=1:kfor i=2:m-1for j=2:m-1w(i,j)=r*u(i-1,j)+r*u(i,j-1)+r*u(i+1,j)+r*u(i,j+1)+(1-4*r)*u(i,j); end end u=w; end toc t=toc umesh(u)交替方向隐式格式（P-R 格式）：11112222,,1,,1,,1,,122111111112222,,1,,1,,1,,122222222n n n n n n n n j l j l j l j l j l j l j l j l n n n n n n n n j l j l j l j l j l j l j l j l u u u u u u u u h hu u u u u u u u h h ττ+++++-+-+++++++++-+-⎧--+-+⎪=+⎪⎪⎪⎨⎪--+-+⎪=+⎪⎪⎩ 将原差分格式化为：1112221,,1,,1,,12111111222,1,,11,,1,/n n n n n n j l j l j l j l j l j l n n n n n n j l j l j l j l j l j lu u u u u u h u u u u u u λλλλλλλτλλλλλλ++++-+-+++++++-+-⎧-+-=++⎪=⎨⎪-+-=++⎩（2+2）（2-2），其中（2+2）（2-2） 121,,112,22,1,2,222222222222222222n nl j nn j l n j m n m lu u u u u u λλλλλλλλλλλλλλλλλλλλλλλλλλλ+++⎛⎫⎪⎛⎫-+--⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪-+-- ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪-+--⎝⎭⎝⎭⎝⎭ ⎪ ⎪⎝⎭-+-⎛-+--+-⎝1211,,1112,22,11,2,222222n n l j n n j l n j m n m lu u uu u u λλλλλλλλλ++++++⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎛⎫⎪ ⎪⎛⎫⎪-⎫⎛⎫ ⎪ ⎪⎪ ⎪ ⎪- ⎪ ⎪⎪ ⎪ ⎪⎪ ⎪⎪ ⎪⎪= ⎪⎪⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪-⎭⎝⎭⎝⎭ ⎪⎪ ⎪⎝⎭⎩代入边界条件，转化为三对角矩阵122,,212,33,1,121,2222222220022222222n nl j n n j l n j m n m lu u u u u u λλλλλλλλλλλλλλλλλλλ++-+-⎛⎫⎪⎛⎫+--⎛⎫⎛⎫⎛⎫⎪ ⎪⎪ ⎪⎪-+-- ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪=+ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪⎪ ⎪ ⎪-+-⎝⎭⎝⎭⎝⎭⎝⎭ ⎪ ⎪⎝⎭+--+1212,,2112,33,11,121,2222220022222n n l j n n j l n j m n m lu uu u u u λλλλλλλλλλλλλ+++++-+-⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎛⎫⎪⎪⎛⎫⎪-⎛⎫⎛⎫⎛⎫ ⎪ ⎪⎪ ⎪ ⎪-- ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎪ ⎪=+ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪-+-⎝⎭⎝⎭⎝⎭⎝⎭ ⎪ ⎪⎝⎭⎩⎪⎪⎪⎪⎪⎪⎪附追赶法源程序：%-------------------------------------------追赶法求解三对角方程组;function x=zg(a,b,c,d)%--------------------------------------------a:方程组系数矩阵A的下对角元素;%--------------------------------------------b:方程组系数矩阵A的主对角元素;%--------------------------------------------c:方程组系数矩阵A的上对角元素;%--------------------------------------------d:追赶法所求方程的右端向量;%--------------------------------------------l:系数矩阵A所分解成的下三角阵L中的下对角元素了l(i); %--------------------------------------------u:系数矩阵A所分解成的下三角阵U中的主对角元素了u(i); n=length(b);u(1)=b(1);y(1)=d(1);for i=1:n-1 %--------------------------追赶法求解之追过程求解Ly=d;l(i)=a(i)/u(i);u(i+1)=b(i+1)-l(i)*c(i);y(i+1)=d(i+1)-l(i)*y(i);endx(n)=y(n)/u(n); %------------------------追赶法求解之赶过程求解Uz=y;for j=n-1:-1:1if u(j)==0break;elsex(j)=(y(j)-c(j)*x(j+1))/u(j);endend%-----------------------------------------------运用P-R差分格式求解二维扩散方程的初边值问题; function pr(ti,h,t)%-------------------------------------------ti:时间步长h:空间步长;k=t/ti+1;m=1/h+1;r=ti/h^2; %------------------------------ r为网格比;w=ones(m,m);u=ones(m,m); %------------------------输入初始值v=ones(m,m);for i=2:m-1for j=2:m-1u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1));endend%------------------------输入用P-R差分格式求解的三对角矩阵b=ones(1,m-2)*(2+2*r);a=-r*ones(1,m-3);c=-r*ones(1,m-3);A=zeros(m-2,m-2); for i=1:m-2A(i,i)=2-2*r; endfor i=1:m-3 A(i,i+1)=r; A(i+1,i)=r; endp=zeros(m-2,1); p(1)=2*r; p(m-2)=2*r; ticfor l=1:kfor i=2:m-1d1=A*u(i,2:m-1)'+p; d1=d1';w(2:m-1,i)=zg(a,b,c,d1); %-------------------------调用追赶法求解 d2=A*w(2:m-1,i)+p;v(i,2:m-1)=zg(a,b,c,d2); %-------------------------调用追赶法求解 end u=v'; end toc t=toc umesh(0:0.1:1,0:0.1:1,u)局部一维格式：11112222,,1,,1,1,,1,22111111112222,,,1,,1,1,,1222211()2222211()222n n n n n n n n j l j lj l j l j l j l j l j l n n n n n n n n j l j l j l j l j l j l j l j l u u u u u u u u h hu u u u u u u u h hττ+++++-+-+++++++++-+-⎧--+-+⎪=+⎪⎪⎪⎨⎪--+-+⎪=+⎪⎪⎩ 将原格式化为：1112221,,1,1,,1,2111111222,1,,1,1,,1(22)(22),/(22)(22)n n n n n n j l j l j l j l j l j l n n n n n n j l j l j l j l j l j l u u u u u u h u u u u u u λλλλλλλτλλλλλλ++++-+-+++++++-+-⎧-++-=+-+⎪=⎨⎪-++-=+-+⎩其中121,1,122,2,1,2,222222222222222222n nl l nn l l n m l n m lu u u u u u λλλλλλλλλλλλλλλλλλλλλλλλλλλ+++⎛⎫⎪⎛⎫-+--⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪-+-- ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪-+--⎝⎭⎝⎭⎝⎭ ⎪ ⎪⎝⎭-+-⎛-+--+-⎝121,1,1112,2,211,2,222222n n j j n n j j n j m n j mu u uu u u λλλλλλλλλ++++++⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎛⎫⎪ ⎪⎛⎫⎪-⎫⎛⎫ ⎪ ⎪⎪ ⎪ ⎪- ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪⎪= ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪-⎭⎝⎭⎝⎭ ⎪⎪ ⎪⎝⎭⎩代入边界条件，转化为三对角矩阵122,1,122,3,1,21,222222220022222222n nl l nn l l n m l n m lu u u u u u λλλλλλλλλλλλλλλλλλλλλλ+++-⎛⎫⎪⎛⎫+--⎛⎫⎛⎫⎛⎫ ⎪ ⎪⎪ ⎪ ⎪-+-- ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪=+ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪⎪ ⎪ ⎪-+-⎝⎭⎝⎭⎝⎭⎝⎭ ⎪ ⎪⎝⎭+--+121,1,2112,3,211,12,2222002222n n j j n n j j n j m n j mu uuu u u λλλλλλλλλλλλλλ+++++-+⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎛⎫⎪⎪⎛⎫⎪-⎛⎫⎛⎫⎛⎫ ⎪ ⎪⎪⎪ ⎪ ⎪-- ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪⎪ ⎪=+ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪-+-⎝⎭⎝⎭⎝⎭⎝⎭ ⎪ ⎪⎝⎭⎩⎪⎪⎪附源程序：%------------------------------------------运用局部一维格式求解二维扩散方程的初边值问题;function god(ti,hi,t) %------------------------------------------------ti 为时间步长 ， hi 为空间步长; m=1/hi;n=t/ti;g=ti/(hi^2); %--------------------------------- g 为网格比 u=ones(m+1,m+1); %------------------------输入初始值 for i=2:m for j=2:mu(i,j)=sin(pi*(i-1)*hi)*sin(2*pi*(j-1)*hi); end enda(1:m-2)=-0.5*g;b(1:m-1)=1+g;c(1:m-2)=-0.5*g; %------------------------输入用局部一维差分格式求解的三对角矩阵 B=zeros(m-1,m+1);for i=1:m-1B(i,i)=0.5*g;B(i,i+1)=1-g;B(i,i+2)=0.5*g;endf=zeros(m-1,1);f(1,1)=0.5*g;f(m-1,1)=0.5*g;w=ones(m+1,m+1);for i=1:nfor j=2:md=B*u(:,j)+f;%-------------------------调用追赶法求解x=zg(a,b,c,d);w(2:m,j)=x';endfor j=2:me=B*w(j,:)'+f;x=zg(a,b,c,e); %-------------------------调用追赶法求解u(j,2:m)=x;endendumesh(u)古典显式在t=1时运行结果：gdxs（0.0025，0.1，1）所用时间t=01.00000 000000 000 1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000 000000 000 0.999999999707990.999999999445260.999999999235500.999999999102410.999999999055020.999999999102400.999999999235500.999999999445260.999999999707991.000000000000001.00000 000000 000 0.999999999445260.999999998943480.999999998547660.999999998290520.999999998204810.999999998290520.999999998547660.999999998943480.999999999445261.000000000000001.00000 000000 000 0.999999999235500.999999998547660.999999997998510.999999997650070.999999997526020.999999997650070.999999997998510.999999998547660.999999999235501.000000000000001.00000 000000 000 0.999999999102400.999999998290520.999999997650070.999999997234010.999999997095320.999999997234010.999999997650070.999999998290520.999999999102401.000000000000001.00000 000000 000 0.999999999055020.999999998204810.999999997526020.999999997095320.999999996941990.999999997095320.999999997526020.999999998204810.999999999055021.00000000000000000000 000 9999102409998290529997650079997234019997095329997234019997650079998290529999102400000000001.00000 000000 000 0.999999999235500.999999998547660.999999997998510.999999997650070.999999997526020.999999997650070.999999997998510.999999998547660.999999999235501.000000000000001.00000 000000 000 0.999999999445260.999999998943480.999999998547660.999999998290520.999999998204810.999999998290520.999999998547660.999999998943480.999999999445261.000000000000001.00000 000000 000 0.999999999707990.999999999445260.999999999235500.999999999102400.999999999055020.999999999102400.999999999235500.999999999445260.999999999707991.000000000000001.00000 000000 000 1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000000000P-R 格式t=1时运行结果：pr（0.0025，0.1，1）所用时间t=0.360000000000001.00000 000000 000 1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000 000000 000 0.999999999544420.999999999133430.999999998807270.999999998597870.999999998525710.999999998597870.999999998807270.999999999133430.999999999544421.000000000000001.00000 000000 000 0.999999999133430.999999998351690.999999997731300.999999997332980.999999997195730.999999997332980.999999997731300.999999998351690.999999999133431.000000000000001.00000 000000 000 0.999999998807270.999999997731300.999999996877400.999999996329170.999999996140260.999999996329170.999999996877400.999999997731300.999999998807271.000000000000001.00000 000000 000 0.999999998597870.999999997332980.999999996329170.999999995684680.999999995462600.999999995684680.999999996329170.999999997332980.999999998597871.000000000000001.00000 000000 000 0.999999998525710.999999997195730.999999996140260.999999995462600.999999995229100.999999995462600.999999996140260.999999997195730.999999998525711.000000000000001.00000 000000 000 0.999999998597870.999999997332980.999999996329170.999999995684680.999999995462600.999999995684680.999999996329170.999999997332980.999999998597871.000000000000001.00000 000000 000 0.999999998807270.999999997731300.999999996877400.999999996329170.999999996140260.999999996329170.999999996877400.999999997731300.999999998807271.000000000000001.00000 000000 000 0.999999999133430.999999998351690.999999997731300.999999997332980.999999997195730.999999997332980.999999997731300.999999998351690.999999999133431.00000000000000000000 000 9999544429999133439998807279998597879998525719998597879998807279999133439999544420000000001.00000 000000 000 1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000000000局部一维格式t=1时的运行结果：god（0.0025，0.1，1）所用时间t= 0.390000000000001.00000 000000 000 1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000 000000 000 0.999999999521570.999999999089960.999999998747440.999999998527530.999999998451750.999999998527530.999999998747440.999999999089960.999999999521571.000000000000001.00000 000000 000 0.999999999089960.999999998269010.999999997617490.999999997199200.999999997055060.999999997199200.999999997617490.999999998269010.999999999089961.000000000000001.00000 000000 000 0.999999998747440.999999997617490.999999996720760.999999996145030.999999995946640.999999996145030.999999996720760.999999997617500.999999998747441.000000000000001.00000 000000 000 0.999999998527530.999999997199200.999999996145030.999999995468210.999999995234990.999999995468210.999999996145030.999999997199200.999999998527531.000000000000001.00000 000000 000 0.999999998451750.999999997055060.999999995946640.999999995234990.999999994989770.999999995234990.999999995946640.999999997055060.999999998451751.000000000000001.00000 000000 000 0.999999998527530.999999997199200.999999996145030.999999995468210.999999995234990.999999995468210.999999996145030.999999997199200.999999998527531.000000000000001.00000 000000 000 0.999999998747440.999999997617490.999999996720760.999999996145030.999999995946640.999999996145030.999999996720760.999999997617500.999999998747441.000000000000001.00000 000000 000 0.999999999089960.999999998269010.999999997617490.999999997199200.999999997055060.999999997199200.999999997617500.999999998269010.999999999089961.000000000000001.00000 000000 000 0.999999999521570.999999999089960.999999998747440.999999998527530.999999998451750.999999998527530.999999998747440.999999999089960.999999999521571.000000000000001.00000 000000 000 1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000000000结论：由上面的表格数据可知：古典显式格式的计算速度最快，且当1/4λ≤时，才是稳定的，局部一维格式的计算速度最慢，但是它是无条件稳定的，P-R格式的速度居中，也是无条件稳定的。

一．上机目的1. 通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；2. 通过上机计算，了解舍入误差所引起的数值不稳定性；3. 熟悉并掌握拉格朗日插值多项式、牛顿插值多项式和分段低次插值，注意其不同特点；4. 了解最小二乘法的基本原理，能通过计算机解决实际问题。

二．上机环境MATLAB 软件等。

三．上机内容1.数值算法稳定性实验；2.插值法实验：拉格朗日插值、牛顿插值以及分段低次插值；3.曲线拟合实验：最小二乘法。

四．实验内容1、数值稳定性实验对n=0,1,2,…20计算定积分dx x x y n n ⎰+=105算法1 利用递推公式151--=n n y ny n=1,2,…,20 取182322.05ln 6ln 51100≈-=+-⎰dx x y 代码：y(1)=log(6)-log(5);for i=1:20y(i+1)=1/i-5*y(i);endk=ones(7,3);for i=1:7for j=1:3k(i,j)=y(3*(i-1)+j);endenddigits(6)vpa(k)结果：[ 0.182322, 0.0883922, 0.0580389][ 0.0431387, 0.0343063, 0.0284684][ 0.0243249, 0.0212326, 0.0188369][ 0.0169265, 0.0153676, 0.0140713][ 0.0129767, 0.0120398, 0.0112295][ 0.0105192, 0.00990388, 0.00930414][ 0.00903483, 0.00745741, 0.012713]算法2 利用递推公式n n y n y 51511-=- n=20,19,.…,1 注意到105151561126110200201020=≤+≤=⎰⎰⎰dx x dx x x dx x 取008730.0)12611051(2120≈+≈y 代码： y(21)=0.008730;for i=2:21j=22-i;y(j)=1/(5*j)-1/5*y(j+1);endk=ones(7,3);for i=1:7for j=1:3k(i,j)=y(3*(i-1)+j);endenddigits(6) ;vpa(k)结果：[ 0.182322, 0.0883922, 0.0580389][ 0.0431387, 0.0343063, 0.0284684][ 0.0243249, 0.0212326, 0.0188369][ 0.0169265, 0.0153676, 0.0140713][ 0.0129766, 0.0120399, 0.0112292][ 0.0105205, 0.0098975, 0.00933601][ 0.00887552, 0.008254, 0.00873]说明：从计算结果可以看出，算法1是不稳定的，而算法2是稳定的。

偏微分方程报告范文偏微分方程是研究多变量函数的微分方程，其中函数的未知量既依赖于自变量，又依赖于多个自变量。

偏微分方程在物理学、工程学和经济学等领域中有着广泛的应用。

本报告将介绍偏微分方程的基本概念、求解方法和应用领域。

一、偏微分方程的基本概念偏微分方程是由未知函数的偏导数和自变量构成的方程。

常见的偏微分方程包括波动方程、传热方程和扩散方程等。

偏微分方程根据阶数可分为一阶和二阶偏微分方程。

一阶偏微分方程中只涉及到未知函数的一阶偏导数，一般可以通过变量分离的方法求解。

二阶偏微分方程中涉及到未知函数的二阶偏导数，求解方法一般包括分离变量法、特征线法和变换法等。

二、偏微分方程的求解方法1.分离变量法：假设未知函数可以表示为两个只依赖于单个自变量的函数的乘积形式，然后将该形式代入到偏微分方程中，再将方程两边关于不同的自变量求积分，从而得到方程的通解。

2.特征线法：通过特征线曲线的方法将偏微分方程转化为常微分方程。

先找出特征线曲线，然后在特征线上引入新的变量，使得偏微分方程变为常微分方程，进而求解。

3.变换法：通过适当的变量变换，将原偏微分方程转化为一个更容易求解的形式。

常用的变换方法有坐标变换、函数变换和变量替换等。

三、偏微分方程的应用领域1.物理学：偏微分方程在物理学中有着广泛的应用。

例如，波动方程可以描述声波、光波和电磁波等在介质中的传播；传热方程可以描述热传导过程；薛定谔方程和波恩-奥本海默方程可以描述量子力学中的粒子行为等。

2.工程学：偏微分方程在工程学中被广泛应用于流体力学、结构力学和电磁场等领域。

例如，纳维-斯托克斯方程用于描述流体的运动；弹性方程用于描述结构的变形和应力分布等。

3.经济学：偏微分方程在经济学中应用较多，尤其是在金融学中。

例如，布莱克-斯科尔斯方程用于定价期权；黑-舒尔斯方程用于描述衍生品的定价和风险管理等。

通过对偏微分方程的研究和求解，可以更好地理解自然界的现象和规律，并为解决实际问题提供数学模型和解决方法。

偏微分方程数值解法上机报告（一）一、实验题目：用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 法观点出发导出的求解处置问题数值解的线性有限元法。

实验一一、实验名称：微分方程数值解二、实验内容一天，缉私舰雷达（A）发现，距离10公里处的海面上有一艘走私船（B）正以匀速25公里/小时沿垂直于AB的直线行驶，缉私舰立即以最大速度（匀速40公里/小时）追赶。

若用雷达进行跟踪，保持船的瞬时速度方向始终指向走私船。

(1) 缉私舰的运动轨迹是怎样的？是否能够追上走私船？如果能追上，需要用多长时间？(2)若缉私舰在走私船2公里处，被走私船发现，走私船立刻加速逃窜（匀速30公里/小时）。

假设走私船沿任意方向逃离，此时是否能够追上走私船？如果能追上，最多需要用多长时间？走私船沿哪个方向逃离时被追上所花的时间最长？(3)如果走私船发现缉私舰后，在逃离的过程中，随时都可能改变方向，此时缉私舰是否能够追上走私船？如果能追上，最多需要用多长时间？三、实验目的熟悉matlab的符号运算环境，掌握微分方程建立、求解方法四、实验原理建立相应的微分方程，并用所掌握的求解方法进行求解与分析.六、实验结果及分析七、实验结论和注记实验二一、实验名称：动态系统模拟二、实验内容一小超级市场有4个付款柜，每个柜台为一位顾客计算货款数的时间与顾客所购商品件数成正比（大约每件费时2s），20%的顾客用支票或信用卡支付，这需要1.5min，其余的顾客付现金则仅需0.5min。

有人倡议设一个快速服务台专为购买8个或8个以下商品的顾客服务，指定另外两个为“现金支付柜”。

请你建立一个模拟模型，用于比较现有系统和倡议的系统的运转。

假设顾客到达平均间隔时间是15秒，顾客购买商品件数按如下频率表分布。

三、实验目的熟悉matlab的统计工具箱，掌握随机数的生成、固定时间增量法与面向事件法的动态模拟方法.四、实验原理把该问题抽象为多队列、多服务台的排队系统，应用动态模拟方法得到模拟数据，分析平均队列长度、平均等待时间等指标。

六、实验结果及分析七、实验结论和注记。

微分方程数值解实验指导（一）实验一: 二阶椭圆型方程差分格式的程序实现1. 实验内容用五点差分格式求解Poisson 方程边值问题⎩⎨⎧∂∈=∈=-=∆,),(,0,),(,:2G y x u G y x f u （1） 其中}1|||,||),{(<=y x y x D 。

（1）用正方形网格)(h k =列出相应的差分方程。

（2）对161,81,51,21=h 分别求出数值解，观察数值解的情况，分析计算结果。

（最好画出数值解的图形）注意：在区域G 的边界上为齐次Dirichlet 条件，在这类边界上不需要给出差分格式。

2. 实验目的及要求按照给定的差分格式编程实现求出数值解；结合格式的相容性和收敛性条件简单分析计算结果。

要求在实验课上算出数值结果；按要求格式写出实验报告；下次实验课前交本次实验的实验报告。

3. 实验原理与实验过程: 以下是求解问题的步骤第一步: 对求解区域作网格剖分。

按照要求得网格剖分图如下（图1为21=h 的情况，图2为51=h 的情况）-1-0.8-0.6-0.4-0.20.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81区域划分示意图图1-1-0.8-0.6-0.4-0.200.20.40.60.81-1-0.8-0.6-0.4-0.200.20.40.60.81区域划分示意图图2第二步: 差分法的目的是：求边值问题的解),(y x u 在节点),(j i y x 的近似值ij u （m j n i ,...,2,1,,2,1==， ）。

为此，需要构造逼近微分方程定解问题的差分格式。

采用五点差分格式，取定x 轴和y 轴方向的步长相等，即h=k ，作两族与坐标轴平行的直线：ih x i +-=1，i=0,1,…,2/h,jh y j +-=1，j=0,1,…,2/h,两族直线的交点为网点（或节点）。

位于G 内的节点为内点，位于G 的边界上的节点为界点。

偏微分方程数值实验报告一实验题目：利用有限差分法求解.0)1(,0)1(),()()(==-=+''-u u x f x u x u 真解为)1()(22x ex 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 段，每个剖分单元的剖分步长记为na b 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 hx 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 .程序代码：第一步：编写有限差分格式相关函数function [ x,U ]=FDld_bvp(N,f,a,b,u)%******************************************************************** %% FD1d_bvp利用中心差分格式求解两点边值问题%参数：% 输入参数：% 整数N,网格节点数，% 函数f(x），计算右端函数f(x);% a，计算区间左端点% b，计算区间右端点% u，真解函数% 输出参数：% 差分解向量U% 均匀剖分区间[a,b]，得到网格x(i)=a+(i-1)*(b-a)/(N-1)h=(b-a)/(N-1);x=(a:h:b)';% 创建线性差分方程组系数矩阵c1=-1/h/h;c2=2/h/h+1;g=[c1*ones(1,N-2),0];c=[0,c1*ones(1,N-2)];d=[1,c2*ones(1,N-2),1];A=diag(g,-1)+diag(d)+diag(c,1);% 创建线性差分方程组右端项rhs=f(x);rhs(1)=u(x(1));rhs(N)=u(x(N));% 求解上述代数系统U=A\rhs;endfunction[e0,e1,emax]=FD1d_error(x,U,u_exact)%% FD1d_ERROR 计算有限差分误差% 参数：% 输入参数：% x，网格节点坐标向量% U，网格节点坐标向量上的有限差分数值解向量Ux% u_exact,真解函数% 输出参数：% e0% e1% emaxN=length(x);h=(x(end)-x(1))/(N-1);ue=u_exact(x);%真解在网格点处的值xee=ue-U;e0=h*sum(ee.^2);e1=sum((ee(2:end)-ee(1:end-1)).^2)/h;e1=e1+e0;e0=sqrt(e0);e1=sqrt(e1);emax=max(abs(ue-U));endfunction FD1d_bvp_test%%测试脚本% 初始化相关数据N=[6,11,21,41,81];L=-1;R=1;emax=zeros(5,1);e0=zeros(5,1);e1=zeros(5,1);%%求解并计算误差for i = 1:5[x,U] =FD1d_bvp(N(i),@f ,L,R,@u);[e0(i),e1(i),emax(i)]=FD1d_error(x,U,@u);X{i}=x;UN{i}=U;endue=u(X{5});%% 显示真阶及不同网格剖分下的数值解plot(X{5},ue,'-k*',X{1},UN{1},'-ro',X{2},...UN{2},'-gs',X{3},UN {3},'-bd ',...X{4},UN{4},'-ch ',X{5} , UN {5},'-mx');title('The solution plot');xlabel('x');ylabel ('u');legend('exact','N=6','N =11','N=21','N =41','N =81'); %% 显示误差format shortedisp ('emax e0 e1 ');disp ([ emax , e0 , e1 ]);end第二步：编写方程的右端函数和真解分别保存为m f .和m u . function f=f(x)f=exp(-x.^2).*(4.*x.^4-15.*x.^2+5);endfunction u=u(x)u=exp(-x.^2).*(1-x.^2);end实验结果：在命令窗口输入>> FD1d_bvp_test回车可得运算结果和图像emax e0 e15.8219e-02 5.3470e-02 1.1724e-011.5919e-02 1.2802e-022.9349e-023.9305e-03 3.1663e-03 7.3357e-039.7959e-04 7.8946e-04 1.8338e-032.4471e-04 1.9723e-04 4.5844e-04。

偏微分中心差分格式实验报告实验目的：1.掌握偏微分的中心差分格式；2.理解中心差分格式的精度和稳定性。

实验原理：中心差分是一种常用的数值求解偏微分方程的格式，其基本思想是用函数在两个点的导数的平均值来近似函数在这两个点中间的导数值。

具体来说，对于一维的偏微分方程，中心差分格式可以表述为：f'(x)=(f(x+h)-f(x-h))/(2h)其中f'(x)表示x点处的导数，h表示步长。

实验步骤：1.编写一个计算函数在任意给定点x处的导数值的中心差分程序；2.给定一个函数f(x)，例如f(x)=x^2，计算在一定范围内的该函数在每个点处的导数值；3.比较计算的导数值与理论值的差异，并分析中心差分格式的精度；4.对给定步长h，逐渐减小h，计算导数值，并观察数值的变化，分析中心差分格式的稳定性。

实验结果与分析：以函数f(x)=x^2为例，给定步长h=0.1，计算在范围[-1,1]内的函数f(x)在每个点处的导数值。

实验结果如下表所示：x，f'(x)，理论值，误差-1.0，-1.999，-2，0.001 -0.9，-1.899，-1.8，0.099 -0.8，-1.698，-1.6，0.098 -0.7，-1.397，-1.4，0.003 -0.6，-0.996，-1，0.004 -0.5，-0.495，-0.5，0.005 -0.4，0.204，0，0.204-0.3，0.615，0.6，0.015 -0.2，1.216，1.2，0.016 -0.1，1.797，1.8，0.003 0.0，1.996，2，0.0040.1，2.193，2.2，0.007 0.2，2.792，2.8，0.008 0.3，3.293，3.3，0.007 0.4，3.594，3.6，0.006 0.5，3.896，3.9，0.004 0.6，4.437，4.4，0.037 0.7，4.998，5，0.0020.9，6.795，6.8，0.0051.0，7.993，8，0.007从实验结果可以看出，随着x的增大，计算的导数值与理论值之间的误差也在增大，但整体上相对较小。

概率论上机实验报告概率论上机实验报告引言：概率论是数学中的一个重要分支，它研究的是随机现象的规律性。

概率论的应用十分广泛，涵盖了自然科学、社会科学、工程技术等各个领域。

为了更好地理解概率论的基本概念和方法，我们进行了一系列的上机实验，通过实际操作来探索概率事件的发生规律以及概率计算的方法。

实验一：硬币抛掷实验在这个实验中，我们使用了一枚标准的硬币，通过抛掷硬币的方式来研究硬币正反面出现的概率。

我们抛掷了100次硬币，并记录了每次抛掷的结果。

通过统计实验结果，我们可以得出硬币正反面出现的频率。

实验结果显示，硬币正面出现的次数为55次，反面出现的次数为45次。

根据频率的定义，我们可以计算出正面出现的概率为55%。

这个结果与我们的预期相符，说明硬币的正反面出现具有一定的随机性。

实验二：骰子掷掷实验在这个实验中，我们使用了一个六面骰子，通过投掷骰子的方式来研究各个面出现的概率。

我们投掷了100次骰子，并记录了每次投掷的结果。

通过统计实验结果，我们可以得出各个面出现的频率。

实验结果显示，骰子的六个面出现的次数分别为15次、18次、17次、16次、19次和15次。

根据频率的定义，我们可以计算出各个面出现的概率分别为15%、18%、17%、16%、19%和15%。

这个结果表明，在足够多次的投掷中，各个面出现的概率是相等的。

实验三：扑克牌抽取实验在这个实验中，我们使用了一副标准的扑克牌，通过抽取扑克牌的方式来研究各个牌面出现的概率。

我们随机抽取了100张扑克牌，并记录了每次抽取的结果。

通过统计实验结果，我们可以得出各个牌面出现的频率。

实验结果显示，各个牌面出现的次数相差不大，都在10次左右。

根据频率的定义，我们可以计算出各个牌面出现的概率都约为10%。

这个结果说明，在足够多次的抽取中，各个牌面出现的概率是相等的。

实验四：随机数生成实验在这个实验中，我们使用了计算机生成的随机数，通过生成随机数的方式来研究随机数的分布规律。

偏微分方程总结报告一、引言偏微分方程是数学中一个重要的分支，它描述了时间和空间中变化的物理量之间的关系。

在自然科学、社会科学和工程学中，偏微分方程有着广泛的应用。

本文将对偏微分方程的基本概念、分类和常见的求解方法进行总结。

二、偏微分方程的基本概念偏微分方程是一个包含未知函数的偏导数的方程。

它通常表示为一个数学表达式，其中包含一个或多个未知函数和这些函数的偏导数。

例如，热传导方程、波动方程和拉普拉斯方程等都是偏微分方程的实例。

三、偏微分方程的分类根据不同的分类标准，偏微分方程可以分为多种类型。

常见的分类方式包括：1. 按照阶数：一阶偏微分方程、二阶偏微分方程等。

2. 按照自变量的个数：常微分方程、偏微分方程等。

3. 按照边界条件：Dirichlet边界条件、Neumann边界条件和Robin边界条件等。

4. 按照方程的形式：线性偏微分方程和非线性偏微分方程等。

四、偏微分方程的求解方法求解偏微分方程的方法有很多种，下面列举几种常见的求解方法：1. 分离变量法：将偏微分方程转化为多个常微分方程，然后求解这些常微分方程。

这种方法适用于具有周期性解的偏微分方程。

2. 有限差分法：将偏微分方程转化为差分方程，然后在离散点上求解这个差分方程。

这种方法适用于具有规则网格的偏微分方程。

3. 有限元法：将偏微分方程转化为变分问题，然后使用有限元方法求解这个变分问题。

这种方法适用于具有复杂边界条件的偏微分方程。

4. 谱方法：将偏微分方程转化为谱问题，然后使用傅里叶分析、小波分析等方法求解这个谱问题。

这种方法适用于具有快速收敛解的偏微分方程。

二阶常微分方程的中心差分求解学校:中国石油大学（华东）理学院 姓名：张道德一、 实验目的1、 构造二阶常微分边值问题:22,(),(),d uLu qu f a x bdx u a u b αβ⎧=-+=<<⎪⎨⎪==⎩其中,q f 为[,]a b 上的连续函数，0;,q αβ≥为给定常数的中心差分格式；2、 根据中心差分格式求解出特定例题的数值解，并与该例题的精确解进行比较。

二、 中心差分格式的构造将区间[,]a b 分成N 等分，分点为： 0,1,2,,i x a ih i N =+=()/h b a N =-。

于是我们得到区间的一个网络剖分。

称为网格的节点称为步长。

得中心差分格式为：11202,1,2,,1,,.i i i h i i i i N u u u L u q u f i N h u u αβ+--+⎧=-+==-⎪⎨⎪==⎩其中式中(),()i i i i q q x f f x ==。

三、 差分格式求解根据中心差分格式可以构造出：1112222222233322212211210012101201001200N N N u f q h h u f q h h h u f q h hh q u f h h ---⎡⎤⎡⎤⎡⎤+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=-+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦可以看出系数矩阵为三对角矩阵，而对于系数矩阵为三对角矩阵的方程组可以用“追赶法”求解，则可以得出二阶常微分方程问题的数值解。

四、 举例求解我们选取的二阶常微分方程边值问题为：222242,01(0)1,(1),x d u Lu x u e x dx u u e ⎧=-+=-<<⎪⎨⎪==⎩其精确解为：2x u e=。

则我们具体求解出的解如下：1、 当10N =时，数值解与精确解为： (1) 表1、10N =时，数值解与精确解统计表x 的值 0.10.20.30.40.5 u 的数值解 1.011069 1.042744 1.096904 1.176896 1.28789 u 的精确解 1.01005 1.040811 1.094174 1.173511 1.284025 两者之差 0.001019 0.001934 0.002729 0.003385 0.003864 x 的值 0.60.70.80.9u 的数值解 1.437443 1.636363 1.900001 2.250209 u 的精确解 1.433329 1.632316 1.896481 2.247908 两者之差0.0041140.0040460.003520.002301将两者绘于同一图中如下：(2)结论：可以看出数值解与精确解之间的误差很小， 在 210-这样一个数量级上。

图像处理偏微分方法实验报告偏微分方程在图像处理领域的研究开始于六、七十年代，从最初的去噪角度和图像恢复的角度出发，相应地引入了偏微分方程。

但直到九十年代才比较系统地将偏微分方程引入图像处理领域，结合其他一些数学工具如数学形态学和仿射几何等，形成了比较完整的理论体系。

由s．Osher等提出的水平集方法具有良好的几何插值性，将灰度插值转化为曲线插值，在图像处理中有很大的影响，在图像修复与去噪、边缘检测、图像匹配、图像识别等方面都取得了相对较好的结果。

一、偏微分方程在图像处理中的应用偏微分方程在图像处理领域的研究开始于六、七十年代，Gabor，Rudin，Osher 等人分别从去噪的角度出发，把偏微分方程引入了图像处理领域。

同时，Osher 还从图像恢复的角度出发，使用最小化全局变差(Total Vadation)方法和变分法，也使用了相应的偏微分方程。

但直到上世纪九十年代，Alvarez，Lions，Morel等才首次系统地将偏微分方程引入了图像处理领域，在理论上有比较大的突破，并结合其他数学工具，如数学形态学、变分法、逼近论、仿射几何等，形成了比较完整的理论体系。

由于理论体系本身的优点，从二维图像(静止图像，如照片)到三维图像(运动图像，如电影)应用的拓展相对容易一些。

这也是偏微分方程在图像处理研究领域中的一个优势。

当然，时间这一维数有它本身的特点，与图像平面有本质的区别，在应用中需要加以注意。

概括地说，偏微分方程的图像处理技术具有以下的特点：(1)基于PDE的方法具有良好的数学基础，可以提供深刻的理论结果，并且算法具备良好的稳定性；(2)一些经典的方法如高斯滤波、中值滤波、膨胀和腐蚀在PDE的统一框架下得到全新的解释；(3)这种视角产生了新的方法，它们比经典的方法包容更多的不变性，如保持结构的滤波、线形增强等；(4)PDE是连续的模型，与具体的离散网格无关，并且具有旋转不变性。

偏微分方程在图像处理中的研究主要集中在图像处理的两个基础部分：图像滤波和图像边缘提取。

2009级数学与应用数学和信息与计算科学专业偏微分方程数值解上机实验实验题目利用有限元方法和有限差分方法求解偏微分方程完成日期 2012年12月17日学生姓名张灵刚所在班级 1102090 任课教师王晓东西北工业大学理学院应用数学系目录一．实验目的 (2)二．实验要求 (2)三．实验题目 (3)四．实验二 (4)1.实验内容 (4)2.实验原理.................................................（4）. 3算法流程. (5)4结果分析 (5)5总结讨论 (6)6源程序 (6)五. 实验三1.实验内容 (17)2.实验原理...............................................（17）. 3算法流程. (18)4结果分析……………………………………….….（18）.5总结讨论 (21)6源程序 (21)偏微分方程数值解上机实验报告实验地点：数学系机房实验时间：第13—15周，周一、四下午5、6节实验分数：占期末考试成绩的30%一、实验目的及意义掌握有限元方法和有限差分方法的程序实现；学会选择合适的有限差分格式求解一维非线性对流占优的非定常对流扩散问题；学会使用三角线性元和四边形线性元的有限元方法求解二维椭圆方程边值问题，并对计算结果进行收敛性分析；尝试采用有限元方法或有限差分方法实现二维初边值抛物型方程的大规模数值求解。

通过实验可以提高学生的动手能力，加深学生对算法的理解。

二、实验要求在下列给出的三个问题中，最少选择两个问题进行编程实现。

要求给出格式的推导过程、算法流程、实现程序、选取的网格参数、以表格或图形的方式给出计算结果、对计算结果进行分析、最后对实验进行总结和讨论。

问题2：用三角线性元和四边形线性元的有限元方法求解方程：28cos(2)sin(2),(,)(0,1)(0,1)(,0)(,1)0;(0,)(1,)sin(2).u x y x y u x u x u y u y y ππππ-∆=∈⨯====取=1/4, 1/8, 1/16, 1/32, 1/64,h 比较两种方法的计算精度，并给出数值收敛率.问题3：选用合适的数值方法求解方程：22122(444),(,)(0,1)(0,1)(0,,)(1,,)(,0,)(,1,)0;(,,0)sin()cos().y y ux y u x y tu y t u y t u x t u x t u x y x y ππππ-∂=-++∆∈⨯∂''=====求0.1,0.5,1.0t =时，点3331(,)6464、1517(,)6464、4749(,)6464、7137(,)128128、4367(,)128128、10989(,)128128、127129(,)256256、6391(,)256256、3133(,)256256处的数值解。

偏微分方程数值解上机实验报告（一）实验一一、上机题目：用线性元求解下列边值问题的数值解：-y′′+π24y=π22sinπ2x，0<x<1y(0)=0，y′(1)=0二、实验程序：function S=bzx=fzero(@zfun,1);[t y]=ode45(@odefun,[0 1],[0 x]);S.t=t;S.y=y;plot(t,y)xlabel('x:´从0一直到1')ylabel('y')title('线性元求解边值问题的数值解')function dy=odefun(x,y)dy=[0 0]';dy(1)=y(2);dy(2)=(pi^2)/4*y(1)-((pi^2)/2)*sin(x*pi/2);function z=zfun(x);[t y]=ode45(@odefun,[0 1],[0 x]);z=y(end)-0;三、实验结果：1．以步长h=0.05进行逐步运算，运行上面matlab程序结果如下：2．在0<x<1区间上，随着x 的不断变化，x ，y 之间关系如下图所示：（二）实验二四、 上机题目：求解Helmholtz 方程的边值问题：21u k u -∆-=，于(0,1)*(0,1)Ω=0u =，于1{0,01}{01,1}x y x y Γ==≤≤≤≤= 12{0,01}{01,1}0,{01,0}{1,01}x y x y u x y x y n Γ==≤≤≤≤=∂=Γ=≤≤==≤≤∂于其中k=1,5,10,15,20五、实验程序：（采用有限元方法，这里对[0,1]*[0,1]采用n*n的划分，n为偶数）n=10;a=zeros(n);f=zeros(n);b=zeros(1,n);U=zeros(n,1);u=zeros(n,1);for i=2:na(i-1,i-1)=pi^2/(12*n)+n;a(i-1,i)= pi^2/(24*n)-n;a(i,i-1)= pi^2/(24*n)-n;for j=1:nif j==i-1a(i,j)=a(i,i-1);else if j==ia(i-1,j-1)=2*a(i-1,i-1);else if j==i+1a(i,j)=a(i,i+1);elsea(i,j)=0;endendendendenda(n,n)=pi^2/(12*n)+n;for i=2:nf(i-1,i)=4/pi*cos((i-1)*pi/2/n)-8*n/(pi^2)*sin(i*pi/2/n)+8*n/(pi^2)*s in((i-1)*pi/2/n);endfor i=1:nf(i,i)=-4/pi*cos(i*pi/2/n)+8*n/(pi^2)*sin(i*pi/2/n)-8*n/(pi^2)*sin((i -1)*pi/2/n);end%b(j)=f(i-1,j)+f(i,j)for i=1:(n-1)b(i)=f(i,i)+f(i,i+1);endb(n)=f(n,n);tic;n=20;can=20;s=zeros(n^2,10);h=1/n;st=1/(2*n^2);A=zeros((n+1)^2,(n+1)^2);syms x y;for k=1:1:2*n^2s(k,1)=k;q=fix(k/(2*n));r=mod(k,(2*n));if (r~=0)r=r;else r=2*n;q=q-1;endif (r<=n)s(k,2)=q*(n+1)+r;s(k,3)=q*(n+1)+r+1;s(k,4)=(q+1)*(n+1)+r+1;s(k,5)=(r-1)*h;s(k,6)=q*h;s(k,7)=r*h;s(k,8)=q*h;s(k,9)=r*h;s(k,10)=(q+1)*h;elses(k,2)=q*(n+1)+r-n;s(k,3)=(q+1)*(n+1)+r-n+1;s(k,4)=(q+1)*(n+1)+r-n;s(k,5)=(r-n-1)*h;s(k,6)=q*h;s(k,7)=(r-n)*h;s(k,8)=(q+1)*h;s(k,9)=(r-n-1)*h;s(k,10)=(q+1)*h;endendd=zeros(3,3);B=zeros((n+1)^2,1);b=zeros(3,1);for k=1:1:2*n^2L(1)=(1/(2*st))*((s(k,7)*s(k,10)-s(k,9)*s(k,8))+(s(k,8)-s(k,10))*x+(s(k,9)-s(k,7))*y);L(2)=(1/(2*st))*((s(k,9)*s(k,6)-s(k,5)*s(k,10))+(s(k,10)-s(k,6))*x+(s (k,5)-s(k,9))*y);L(3)=(1/(2*st))*((s(k,5)*s(k,8)-s(k,7)*s(k,6))+(s(k,6)-s(k,8))*x+(s(k ,7)-s(k,5))*y);for i=1:1:3for j=i:3d(i,j)=int(int(((((diff(L(i),x))*(diff(L(j),x)))+((diff(L(i),y))*(dif f(L(j),y))))-((can^2)*L(i)*L(j))),x,0,1),y,0,1);d(j,i)=d(i,j);endendfor i=1:1:3for j=1:1:3A(s(k,(i+1)),s(k,(j+1)))=A(s(k,(i+1)),s(k,(j+1)))+d(i,j);endendfor i=1:1:3b(i)=int(int((L(i)),x,0,1),y,0,1);B(s(k,(i+1)),1)=B(s(k,(i+1)),1)+b(i);endendM=zeros((n+1)^2,n^2);j=n^2;for i=(n^2+n):-1:1if ((mod(i,(n+1)))~=1)M(:,j)=A(:,i);j=j-1;else continueendendpreanswer=M\B;answer=zeros((n+1)^2,1);j=1;for i=1:1:(n^2+n)if ((mod(i,(n+1)))~=1)answer(i)=preanswer(j);j=j+1;else answer(i)=0;endendZ=zeros((n+1),(n+1));for i=1:1:(n+1)^2s=fix(i/(n+1))+1;r=mod(i,(n+1));if(r==0)r=n+1;s=s-1;elseendZ(r,s)=answer(i);end[X,Y]=meshgrid(1:-h:0,0:h:1);surf(X,Y,Z);toc;t=toc;K=a ;B=b';U=inv(K)*Bfor i=1:nu(i,1)=4/(pi^2)*sin(pi*i/n/2);endue=U-u六、实验结果：程序中的变量can为问题中的k，为了方便比较，采用了画图的方式。