偏微分方程解的几道算例(差分、有限元)+含matlab程序

2.时间、空间均为 0 − 2π ，且网格为 20 × 20 的图像结果（数据太多--略去）：
>>u=PP(0,2*pi,0,2*pi,2*pi/20,2*pi/20); >>x=[0:2*pi/20:2*2pi]; >>y=[0:2*pi/20:2*2pi]; >>mesh(x,y,u)
（三）结果分析：
程序 2 ： function u=PP(t0,tn,x0,xn,ht,hx) % t0时间起点 %tn时间终点 %x0空间起端点 %xn空间终端点 %ht为时间步长 nt=(tn-t0)/ht; nx=(xn-x0)/hx; u=zeros(nt+1,nx+1); r=ht/hx; for i=1:(nx+1)
（三）程序（附最后）
实验内容 3：
用线性元求解下列问题的数值解：
⎧∆u = −2,
−1 < x, y <1,
⎪⎨u(x, −1) = u (x ,1) = 0,
−1 < x <1,
⎪⎩ux (−1, y ) = 1,ux (1, y ) = 0,
−1 < x < 1.
（一）算法描述：
（二）实验结果：
A(i,i)=r2; if i>1
A(i-1,i)=-r; A(i,i-1)=-r; end end u=zeros(N+1,M+1); u(N+1,:)=u1; for k=1:N b=u(N+2-k,2:M)+0.02; u(N+1-k,2:M)=inv(A)*b';%求解迭代方程组 end uT=u(1,:);%0.25时刻的解 %精确解与数值解画图 x1=[0,x,1]; plot(x1,uT,'o') hold u_xt = exp (-pi*pi*T)*sin (pi*x1) + x1.*(1 - x1); plot (x1, u_xt, ' r') e=u_xt-uT; 六点格式 function [e]=six(dx,dt,T) %用六点对称格式求解，dx为x方向步长，dt为t方向步长 % e为误差 M=1/dx; N=T/dt; %得到第一层的值 u1=zeros(1,M+1); x=[1:M-1]*dx; u1([2:M])= sin(pi*x)+x.*(1 - x); %网比 r=dt/dx/dx;r2=2+2*r;r3=2-2*r; %构造三对角矩阵A for i=1:M-1 A(i,i)=r2;
a(i*m+1,i*m+1)=4; a((i+1)*m,(i+1)*m)=4; end clear i for j=2:(m-1) a(j,j)=4;
for i=1:m x(i+(k-1)*m)=xl+(k-1)*(xr-xl)/(n-1); y(i+(k-1)*m)=yd+(i-1)*(yu-yd)/(m-1);
end end clear k i a(m,m)=2; a((n-1)*m+1,(n-1)*m+1)=2; a(1,1)=3; a(m*n,m*n)=3; for i=1:(n-2)
《偏微分方程数值解》 上机报告
实验内容 1：
分别用向前差分格式、向后差分格式及六点对称格式，求解下列问题：
⎧ ∂u ∂2u
⎪⎪ ∂t
= ∂x 2
+ 2,
⎨⎪u(0,t) = u(1,t) = 0,
0 < x < 1,t > 0, t > 1,
⎪⎩u(x, 0) = sin(π x ) + x (1− x ).
b(i+1)=r*a(i+2)+(1-2*r)*a(i+1)+r*a(i)+2*t; End 向后格式 function [e]=back(dx,dt,T) %用向后差分格式求解，dx为x方向步长，dt为t方向步长 % e为误差 M=1/dx;
N=T/dt; %得到第一层的值 u1=zeros(1,M+1); x=[1:M-1]*dx; u1([2:M])= sin(pi*x)+x.*(1 - x); %网比 r=dt/dx/dx;r2=1+2*r; %构造三对角矩阵 for i=1:M-1
1.区域 [−1,1]×[−1,1]被剖分成 10×10时的数值和图像结果
>>FE(-1,1,-1,1,10,10) ans=
（结果采用图片截得）
（相应的网格剖分情况） >>u=FE(-1,1,-1,1,10,10);
>>x=[-1:0.2:1]; >>y=[-1:0.2:1]; >>mesh(x,y,u)
x 方向 h = 0.1 ， t方向τ = 0.01.在 t = 0.25 时观察数值解与精确解
u = e−π 2 sin(π x) + x(1− x) 的误差. （一）算法描述：
（二）实验结果：
1.误差的数值解结果数值对比
(A)“向前差分格式”程序：
>>forward(0.1, 0.01, 0.25) Current plot held ans = 0.0000 0.0027 0.0051 0.0082 0.0070 0.0051 (B)“向后差分格式”程序： >>back(0.1, 0.01, 0.25) Current plot held ans = 0.0000 -0.0037 -0.0071 - 0.0114 -0.0097 -0.0071 (C)“六点差分格式”程序： >>six(0.1, 0.01, 0.25) Current plot held ans = 0.0000 -0.0005 -0.0009 -0.0015 - 0.0013 -0.0009
u = for_climb (u, r, t); End %精确解与数值解画图 plot (x, u, 'o') hold u_xt = exp (-pi*pi*T)*sin (pi*x) + x.*(1 - x); plot (x, u_xt, ' r') e=u_xt-u; function b=for_climb(a,r,t) %向前差分格式中，由下一层得到上一层 L=length(a); b=zeros(1,L); for i=1:L-2
本人精通MATLAB等编程语言，可以提供以下方向的帮助 1. MATLAB/GUI/SIMULINK/C++/VC++编程问题； 2. 线性与非线性控制、智能控制、模糊控制； 3. 数值计算问题、小波分析算法、有限元问题； 4. 电机控制、电力系统、机器人路径优化、机器人控制； 5. 粒子群算法、神经网络、模拟退火算法等智能优化算法； 6. 图像处理、信号处理、语音信号处理、电子通信等方向；
程序 3： function w=FE(xl,xr,yd,yu,M,N) %x1为区域左边界值（本题为-1） %xr为区域右边界值（本题为1） %yd为区域下边界值（本题为-1） %yu为区域下边界值（本题为-1） %M为y轴方向剖分格数 %N为X轴方向剖分格数 %网格剖分，单元编号 m=M+1;n=N+1; a=zeros(m*n);b=zeros(m*n); x=zeros(1,m*n);y=zeros(1,m*n); for k=1:n
u
0 j
=
0
，
j
=
0,1,....nx
,
∂u ∂t
|(0, j ) =
u
−1 j
−
u
0 j
ht
= sin x j ,
即u−j 1＝u0j－ht sin xj ,
将上式代入到差分格式可以求得 u1j = ht sin x j ， j = 0,1,.....nx ,
最后在迭代式中利用
u
0 j
，
u1j
，可以求得
有问题的朋友，可以将问题直接发到我的邮箱，24小时内给您答复！非 常欢迎大家加我为QQ好友，欢迎访问我的空间！
联系方式： QQ：626815632 邮箱：626815632@ QQ空间：/
声明：本资料来源于网络，切勿用做商业用途！请您支持正版图书！
注：曲线表示精确解，"o"表示数值解（t=0.25 时）. “六点差分格式”：
注：曲线表示精确解，"O"表示数值解（t=0.25 时）. (三)结果分析
通过（一），（二），我们检验了三种方法都能很好的求解此一维热传导方程，其 中明显能发现“六点对称格式”的误差更小。 （四）程序（附最后）
实验内容 2：
if i>1 A(i-1,i)=-r; A(i,i-1)=-r;
end end %构造三对角矩阵B,方便得到迭代方程组的右端b for i=1:M-1
B(i,i)=r3; if i>1
B(i-1,i)=r; B(i,i-1)=r; end end u=zeros(N+1,M+1); u(N+1,:)=u1; for k=1:N b=B*(u(N+2-k,2:M))'+0.04; u(N+1-k,2:M)=inv(A)*b;%求解迭代方程组 end uT=u(1,:);%0.25时刻的解 %精确解与数值解画图 x1=[0,x,1]; plot(x1,uT,'o') hold u_xt = exp (-pi*pi*T)*sin (pi*x1) + x1.*(1 - x1); plot (x1, u_xt, ' r') e=u_xt-uT;