[关闭]
@joyphys 2015-08-13T06:57:00.000000Z 字数 1811 阅读 3008

Matlab数值计算最简单的一个例子——指数衰减

Blog


放射性衰变是指数衰减的典型例子。另外还有化学反应某反应物的减少,RC电路电流的减小,大气压随海拔高度的减小等。

指数衰减的方程:

dN(t)dt=N(t)τ(1)

其中,N(t)t时刻的物理量N,对于放射性衰变,N就是未衰变的原子核数目。τ为时间常数。

方程(1)有解析解:

N(t)=N(0)exp(t/τ)

这个解可以通过Matlab符号计算求得:

  1. dsolve('DN=-N/tau')
  2. ans =
  3. C3*exp(-t/tau)

数值求解方程(1),可用欧拉格式将方程离散化。

ti=(i1)Δt,i=1,2,,npoints

dN(t)dtN(t)N(tΔt)Δt

将以上两式带入方程(1),得离散之后的方程:
N(ti+1)=N(ti)N(ti)Δtτ

代入初始条件,即可得解。

下面写个Matlab 脚本文件,重复出Computational Physics_Giordano 2nd Edition的图1.1,pp11

  1. %
  2. % Exponent decay
  3. % 'Computational Physics' book by N Giordano and H Nakanishi
  4. % Section 1.2 p2
  5. % Solve the Equation dN/dt = -N/tau
  6. % by Joyful Physics Blog
  7. % ------------------------------------------------------------
  8. N_nuclei_initial = 100; %initial number of nuclei
  9. npoints = 101; % Discretize time into 100 intervals
  10. dt = 0.05; % set time step
  11. tau=1; % set time constant
  12. N_nuclei = zeros(npoints,1); % initializes N_nuclei, a vector of dimension npoints X 1,to being all zeros
  13. time = zeros(npoints,1); % this initializes the vector time to being all zeros
  14. N_nuclei(1) = N_nuclei_initial; % the initial condition, first entry in the vector N_nuclei is N_nuclei_initial
  15. time(1) = 0; %Initialise time
  16. for step=1:npoints-1 % loop over the timesteps and calculate the numerical solution
  17. N_nuclei(step+1) = N_nuclei(step) - (N_nuclei(step)/tau)*dt;
  18. time(step+1) = time(step) + dt;
  19. end
  20. % calculate analytical solution below
  21. t=0:0.05:5;
  22. N_analytical=N_nuclei_initial*exp(-t/tau);
  23. % Plot both numerical and analytical solution
  24. plot(time,N_nuclei,'ro',t,N_analytical,'b'); %plots the numerical solution in red and the analytical solution in blue
  25. xlabel('Time (s)')
  26. ylabel('Number of nuclei')
  27. text(2,80,'Time constant = 1s')
  28. text(2,70,'Time step = 0.05s')

运行程序,得到:

添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注