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! :)

[Sean de Wolski]

There is no z in the above equation.

[Teja Muppirala]

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).

Answer 1

[Teja Muppirala]

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.

[Abhinand Jha]

What is the logic behind creating a First derivative operator?

By which mathematical way did you create it?

Answer 2

[Abhinand Jha]

C''= double derivative w.r.t. z and ofcourse the boundary condition is dc/dz=0 (dc/dx was a typo) and thanks teja