24小时热门版块排行榜

返回列表

zwb565055403

新虫 (小有名气)

应助: 0 (幼儿园)
金币: 86.3
散金: 40
红花: 1
帖子: 88
在线: 112.5小时
虫号: 2783017
注册: 2013-11-06
专业: 偏微分方程

[交流] 编程为什么结果出来不对啊

我写的程序如下（虽然能运行出来，结果却不对大家帮我看下程序是否有问题）：
% Example 1
clear;clc
%%%%%%%%%%%%%%%%%%%% 方程里的参数设置 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha=0.5;
beta=0.5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L=1; %区间x的长度[0,1]
Mx=L/0.0025;  %区间x的等分的份数
h=L/Mx;  %区间x的划分小区间长度
x=[0:Mx]*h;  % x的点值
%上面给出了x的点值，结合t点的点值和函数 u(x,t)的值就可以画出曲面u(x,t),三维数组
N=1000; %时间t的等分份数
tau=0.0025; % 时间t的划分小区间长度
t=[0:N]*tau;  %时间t的点值 0点和1点 101个
  %%%%%%%%%%%%%%%%%%%%%%%%%%% 系数矩阵里面的常数定义 r R(i) %%%%%%%%%%%%%%%%%%
R=(tau^alpha)*(gamma(2-alpha))/(h*h); % r常数
%计算 r(j)
for j=1:Mx-1
  r(j)=beta*R/j;  % 变数 r_j
end
r;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%分数阶系数矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v0=1+2*R-r;  %主对角线上元素
A=diag(v0);
v1=r-R ;
v1=v1(2:end);%下次对角线上的元素
B=diag(v1,-1);
v3=-R*ones(1,Mx-2); %上次对角线上的元素
C=diag(v3,1);
D=A+B+C;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%非齐次项f(x,t)=2\sqrt(t/pi) sin(pi x)-pi t ((cos pi x)/(2x)-pi sin(pi x))
for k=1:size(t,2)    %size(a,2)指行向量元素个数
F(k,

=2*(t(k)^(1/2)/pi^(1/2)).*sin(pi*x) + pi*t(k).*(pi*sin(pi.*x) - cos(pi*x)./(2*x));
%第k行向量时间不变而空间位置变化
end
F=F(2:end-1,2:end-1);
F=F';
%去掉第一行和最后一行去掉第一列和最后一列对应 t=0 x=0
%%%%%%%%%%%%%%%%%%下面用递推关系式得到关于 u 的数表 %%%%%%%%%%%%%%%%%%%%%%%%%%%
u(:,1)=inv(D)*((tau^alpha)*gamma(2-alpha)*F(:,1));
% D*u1=u0+ tau^alpha)F1 例一中 u0=0 即初值函数为0函数
%定义 w(k)=(1+k)^(1-alpha)-k^(1-alpha);d(k)=w(k)-w(k+1)
for k=1:N
w(k)=(1+k)^(1-tau)-k^(1-tau);
end
for k=1:N-1
d(k)=w(k)-w(k+1);    %系数 d(k)
end
d;
%%%%%%%%%%%%%%%%%%%%%%%%%% 单独算出 u2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u(:,2)=inv(D)*((1-w(1))*u(:,1)+(tau^alpha)*gamma(2-alpha)*F(:,1));
%%%%%求和 \sum d(k)*u^(n-k) 列和 sum(A(:,a:b);2)
for n=2:N-2
for k=1:n-1
      s(:,k)=d(n-k)*u(:,k);
end
u(:,n+1)=inv(D)*((tau^alpha)*gamma(2-alpha)*F(:,n+1)+(1-w(1))*u(:,n)+sum(s(:,1:n-1),2));
% u0=0
end
h=x(2:end-1);
g=t(2:end-1);
[h,g]=meshgrid(h,g);%生成网格
mesh(h,g,u');%画出曲面图
xlabel('x')  %添加坐标轴
ylabel('t')
zlabel('u(x,t)')
for k=1:N-1
for j=1:Mx-1
plot3(x(j),t(k),u(j,k))
e(j,k)=u(j,k)-t(k)*sin(pi*x(j));  %误差函数
end
end
e=max(max(abs(e)))    %最大误差

回复此楼