finite differences scheme
显示 更早的评论
Hello,
I have a 1-D steady state (dc/dt=0) differential equation in the atmosphere. It looks like follows,
K*C'' + (K'+K/H)*C' + (1/H*K'- (K/H^2)*H'- (L+Si))C + S = 0
where, C = concentration of the contaminant in the atmosphere at different heights z K = vertical diffusion coefficient H = scale height L = decay constant Si= constant S = source term
C'' = double derivative of C w.r.t. z C',K',H'= derivative of C,K,H w.r.t. z
K,H,Si,S,L are all known values.
I am trying to solve the above differential equation numerically by means of finite differences of 1st order with boundary conditions,
At the top boundary: C = S/L at the bottom boundary k*dc/dx = 0
Can anyone tell me how to write this routine in matlab?
help would be greatly appricieted! :)
2 个评论
Sean de Wolski
2011-5-11
There is no z in the above equation.
Teja Muppirala
2011-5-11
There doesn't need to be a z in that equation. z is the independent variable. C = C(z).
(I think the boundary condition is supposed to be dC/dz = 0).
回答（2 个）
Teja Muppirala
2011-5-11
Here is an example:
Solve
t*y' + (1+sin(t))*y = cos(t)
y(0) = 1, using first order finite differences
N = 2000;
t = linspace(0,1,N)';
dt = t(2) - t(1);
% Create the first derivative operator
D1 = diag(ones(N,1));
D2 = diag(-ones(N,1),-1);
D2 = D2(1:N,1:N);
D = (D1+D2)/dt;
% Write it as M*y = f
M = [diag(t)*D + diag(1+sin(10*t))];
f = cos(t);
% Add boundary conditions
f(1) = 1;
M(1,:) = 0;
M(1,1) = 1;
%Solve and plot
y = M\f;
plot(t,y,'r');
hold on;
% Compare it with MATLAB'S ODE45 solver
[t2,y2] = ode45(@(t,x) (cos(t) - (1+sin(10*t)).*x)./t, [eps 1], 1);
plot(t2,y2,'b:');
legend({'Finite difference solution', 'ODE45 solution'});
Now your problem is a second order differential equation, and what I called y and t, you are calling C and z, but the process is exactly the same. You just need to find a matrix to represent the 2nd derivative operation. I'll leave that part to you.
类别
在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
2011-5-11
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
