steva
2008-01-31, 12:00
请教各位高手.
我在程序里设定了tmax为1E6,但是程序运行结束后,结果显示tmax小了一个数量级,为1E5了,图像的横坐标也是1E5,这是为什么呢?另外我想把解析解和数值解都plot出来,可是解析解为什么出不来呢?恳请各位高手帮我看看这个程序错误在哪个地方。先谢了阿
syms P0 P a b c d y E lamda N t L f n tmax k1 k2 delta_t C1 N_exact
a=input('Input the parameter for production rate of cosmogenic nuclides a:');
b=input('Input the parameter for production rate of cosmogenic nuclides b:');
c=input('Input the parameter for production rate of cosmogenic nuclides c:');
d=input('Input the parameter for production rate of cosmogenic nuclides d:');
y=input('Input the altitude of the sample y:');
E=input('Input the erosion rate E:');
lamda=input('Input the decay constant lamda:');
L=input('Input the production depth of cosmogenic nuclides L:');
delta_t=input('Input the value of the stepsize delta_t: ');
P=a+b*y+c*y^2+d*y^3;
C1=E/L+lamda;
tmax=1E6;
t=0:delta_t:tmax;
f=@(t,N) P-C1*N;
n=tmax/delta_t;
t=0; N=0;
for i=1:n
k1=delta_t*f(t(i),N(i));
k2=delta_t*f(t(i)+delta_t,N(i)+k1);
N(i+1)=N(i)+k1/2+k2/2;
end
N_exact=P/C1-exp(-C1*t)*P/C1;
plot(t,N,'r',t,N_exact,'b')
我在程序里设定了tmax为1E6,但是程序运行结束后,结果显示tmax小了一个数量级,为1E5了,图像的横坐标也是1E5,这是为什么呢?另外我想把解析解和数值解都plot出来,可是解析解为什么出不来呢?恳请各位高手帮我看看这个程序错误在哪个地方。先谢了阿
syms P0 P a b c d y E lamda N t L f n tmax k1 k2 delta_t C1 N_exact
a=input('Input the parameter for production rate of cosmogenic nuclides a:');
b=input('Input the parameter for production rate of cosmogenic nuclides b:');
c=input('Input the parameter for production rate of cosmogenic nuclides c:');
d=input('Input the parameter for production rate of cosmogenic nuclides d:');
y=input('Input the altitude of the sample y:');
E=input('Input the erosion rate E:');
lamda=input('Input the decay constant lamda:');
L=input('Input the production depth of cosmogenic nuclides L:');
delta_t=input('Input the value of the stepsize delta_t: ');
P=a+b*y+c*y^2+d*y^3;
C1=E/L+lamda;
tmax=1E6;
t=0:delta_t:tmax;
f=@(t,N) P-C1*N;
n=tmax/delta_t;
t=0; N=0;
for i=1:n
k1=delta_t*f(t(i),N(i));
k2=delta_t*f(t(i)+delta_t,N(i)+k1);
N(i+1)=N(i)+k1/2+k2/2;
end
N_exact=P/C1-exp(-C1*t)*P/C1;
plot(t,N,'r',t,N_exact,'b')