|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
function My_Experiment_2()
clc;
clear;
b=1;
L=pi;
m=20;
n=1000;
dt=4e-4;
dx=L/m;
alpha=0.5;
d=(2*(1:n+1).^(1-alpha)-(2:n+2).^(1-alpha)-(0:n).^(1-alpha));
Nmax=100;
tol=1e-6;
u=zeros(m+1,n+1);
r=b^2*dt^alpha /(dx^2)*gamma(2-alpha);
c=(1:n+1).^(1-alpha)-(0:n).^(1-alpha);
u(1,1:(n+1))=0;
u(m+1,1:(n+1))=0;
u(2:m,1)=feval('myfun1',dx:dx:(m-1)*dx)';
u1=u;
for k=1:Nmax
for j=3:(n+1)
for i=2:m
u1(i,j)=(c(j)*u(i,1)+d(1:j-1)*u(i,j-1:-1:1)'+r*u(i+1,j)+r*u(i-1,j))/(1+2*r);
end
end
if(round(k/10)==k/10)
fprintf('k=%4.2d\n',k);
end
if(max(abs(u1-u))/max(abs(u1))<tol)
u=u1;fprintf('迭代次数为 :%4.1d',k);
mesh(0:dt:n*dt,0:dx:m*dx,u);
plot(0:dx:m*dx,u(:,end)');
t=0.4;
sum=1;
for i=1:100
sum=sum+(-1).^i*t^(i*alpha)/gamma(i*alpha+1);
end
q=sin(0:dx:m*dx).*sum;
hold on;
plot(0:dx:m*dx,q,'*');
hold off;
return;
end
u=u1;
end
u(5,5)
刚接触MATLAB,好多都不懂啊,请问
u=zeros(m+1,n+1); 这个语句的作用是什么?
改成u=ones(m+1,n+1); 为什么输出的值是不一样的
[ 本帖最后由 eight 于 2007-12-10 21:20 编辑 ] |
|