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

Matlab数值计算最简单的⼀个例⼦——指数衰减放射性衰变是的典型例⼦。

另外还有化学反应某反应物的减少，RC电路电流的减⼩，⼤⽓压随海拔⾼度的减⼩等。

指数衰减的⽅程：
dN(t)
dt=−N(t)τ
其中，N(t)为t时刻的物理量N，对于放射性衰变，N就是未衰变的原⼦核数⽬。

τ为时间常数。

⽅程(1)有解析解：
N(t)=N(0)exp(−t/τ)
这个解可以通过Matlab符号计算求得：
dsolve('DN=-N/tau')
ans =
C3*exp(-t/tau)
数值求解⽅程(1)，可⽤将⽅程离散化。

t i=(i−1)Δt,i=1,2,…,npoints
dN(t) dt≈N(t)−N(t−Δt)
Δt
将以上两式带⼊⽅程(1)，得离散之后的⽅程：
N(t i+1)=N(t i)−N(t i)Δt τ
代⼊初始条件，即可得解。

下⾯写个Matlab 脚本⽂件，重复出Computational Physics_Giordano 2nd Edition的图1.1，pp11 %
% Exponent decay
% 'Computational Physics' book by N Giordano and H Nakanishi
% Section 1.2 p2
% Solve the Equation dN/dt = -N/tau
% by Joyful Physics Blog
% ------------------------------------------------------------
N_nuclei_initial = 100; %initial number of nuclei
npoints = 101; % Discretize time into 100 intervals
dt = 0.05; % set time step
tau=1; % set time constant
N_nuclei = zeros(npoints,1); % initializes N_nuclei, a vector of dimension npoints X 1,to being all zeros
time = zeros(npoints,1); % this initializes the vector time to being all zeros
N_nuclei(1) = N_nuclei_initial; % the initial condition, first entry in the vector N_nuclei is N_nuclei_initial
time(1) = 0; %Initialise time
for step=1:npoints-1 % loop over the timesteps and calculate the numerical solution
N_nuclei(step+1) = N_nuclei(step) - (N_nuclei(step)/tau)*dt;
time(step+1) = time(step) + dt;
end
% calculate analytical solution below
t=0:0.05:5;
N_analytical=N_nuclei_initial*exp(-t/tau);
% Plot both numerical and analytical solution
plot(time,N_nuclei,'ro',t,N_analytical,'b'); %plots the numerical solution in red and the analytical solution in blue xlabel('Time (s)')
ylabel('Number of nuclei')
text(2,80,'Time constant = 1s')
text(2,70,'Time step = 0.05s')
运⾏程序，得到：
Processing math: 100%。