Finite difference method - Second order equation with special conditions

Hello I am trying to solve this problem with the finite difference method.

(1)

Where and

with the following conditions:

(2)

(boundary condition when r = 2)

and this discretization when r = 2 (and only then!)

(3)

My attempt so far is to use central difference approximation with (1) and using (3) instead of first order derivative only for the last approximation (N=Nend). But when I write down the matrix I only get 0's in the b-vector (Ay = b) which feels very wrong. I have no clue how I should use (2) either. Can someone help me to think in the right direction here? Thank you.

[Torsten]

You are not allowed to use bvp4c ?

What is your boundary condition at r=1 ?

[Casper Hae]

I don't know what bvp4c is and we've never used it. But it says in the task that we should use finite difference method

The boundary condition at r=1 is

Answer 1

[Torsten]

As your equation can be written as

1/r * d/dr (r*dT/dr) = 0

the usual discretization in r = r_i is

1/r_i * [(r_i + dr/2)*(T(r_(i+1))-T(r_i))/dr - (r_i - dr/2)*(T(r_i)-T(r_(i-1)))/dr] = 0,

thus

(ri + dr/2)*(T(r_(i+1))-T(r_i)) - (r_i - dr/2)*(T(r_i)-T(r_(i-1))) = 0 (2 <=i <= n-1)

The boundary condition at r = 1 gives

T(r_1) = T_a

The discretization at r_n of the differential equation gives

r_n*dT/dr(r_n) - (r_n-dr/2)*dT/dr(r_n-dr/2) = 0

Inserting

dT/dr(r_n) = -alpha/k*(T(r_n)-T_b)

from the boundary condition and discretizing

dT/dr(r_n-dr/2) = (T(r_n)-T(r_(n-1))/dr

leads to

r_n * alpha/k*(T_b-T(r_n)) - (r_n - dr/2)*(T(r_n)-T(r_(n-1))/dr = 0

or

T(r_n)*r_n*alpha/k + (r_n - dr/2)*(T(r_n)-T(r_(n-1))/dr = r_n*alpha/k*T_b

Thus the vector of the right hand side is not zero - it has entries at the first and last position coming from the boundary values.

[Casper Hae]

Why do you take (1/r) * r when the equation is r? And I don't really understand how you discreticize. Isn't it possible to use central approximation for the second and first derivative? Because they both have an error of O(h^2).

[Torsten]

Why do you take (1/r) * r when the equation is r?

I don't understand.

And I don't really understand how you discreticize. Isn't it possible to use central approximation for the second and first derivative? Because they both have an error of O(h^2).

There are many ways how to discretize the equation. Take the one you like most - there is no right or wrong.

I like the method above because you only need to discretize one expression

( 1/r * d/dr (r * dT/dr ) )

instead of two

( d^2T/dr^2 + 1/r * dT/dr )

If you want to use your method, then only take from my answer: The vector of the right-hand side is not zero, but has entries in the first and last position because of the boundary values.

[Casper Hae]

Yeah I tried my method and that was "approximately" (don't know if I made everything correctly though) what I got as well. So if T(r) represents temperature, will the temperature be 0 everywhere except on the boundaries? Because that seems very weird in the problem given.

[Torsten]

So if T(r) represents temperature, will the temperature be 0 everywhere except on the boundaries?

No.

The general solution is

T( r) = a*log(r ) + b

You can determine a and b from the boundary conditions.

[Casper Hae]

How did you find the general solution? I don't know if I have missed something important.

[Torsten]

   1/r * d/dr (r*dT/dr) = 0
-> d/dr (r * dT/dr) = 0
-> r * dT/dr = a = const
-> dT/dr = a/r 
-> T = a*log(r ) + b

[Casper Hae]

Are you able to show the approximations for central differences? I don't really understand you approximations. I guess I have to use Neumann boundary conditions for the outer boundary, but I don't know how.

[Torsten]

The approximation of the differential equation in inner grid points is done by central differences.

The approximation of the differential equation in r = R is done by backward Euler for the second derivatives, followed by central differencing for the resulting first-derivative term. The dT/dr at r=R is then substituted by the boundary condition.

What exactly is unclear ? I think I cannot write out what to do more clearly than I have done above.

Do you have to use a different discretization ? Then show your attempt.

[Casper Hae]

This is for the inner points, right? (Just an example with different N's).

Then using the boundary conditions at eq. 7 we get the first element i b-vector is To/2h - To*1/h (r=1).

Then the last boundary condition (eq 10.) is very complicated in my calculations. But does it seem true so far?

[Torsten]

Looks fine to me.

But there is no need to solve the first equation for T0.

Just add the equation

T0 = Ta

to your linear system of equations you have to solve.

For the equation in the last grid point, proceed similiarily:

Your last equation reads

r_n/h^2 (T_(n+1) - 2*T_n + T_(n-1)) + 1/(2*h) * (T_(n+1) - T_(n-1)) = 0

and from the boundary condition you get

-k* (T_(n+1)-T_(n-1))/(2*h) = alpha*(T_n - Tb)

Just leave both equations as they are and include them in your linear system of equations.

T_(n+1) is the temperature in the "fictuous" grid point 2+h. (One usually assumes that the heat equation continues to be valid also at the boundary of the domain).

This will give you a system of (n+2) equations for the (n+2) unknowns T0,T1,...,T_(n+1) (I know, T0 is not unknown, but althoug it's known it's simpler to include it in the vector of unknowns as you saw above).

[Casper Hae]

[Casper Hae]

[Casper Hae]

I don't know. Is this right? Last boundary gets quite complicated so I don't know if there's any errors since I am not used to do Neumann boundary conditions.

[Torsten]

Why don't you follow my advice and write out the system for T0,T1,...,T_(n+1) instead of T1,...,T_n ?

The system will be much more clearly structured.

The first two equations from sheet 1 (with n replaced by i in the first equation) are sufficient to deduce the matrix and the vector of the right-hand side as follows:

row 1 of the matrix:       1 0 0 0 0 0 ... 0;...
row 2 of the matrix:       r_1/h^2 - 1/(2*h) , -2*r_1/h^2, r_1/h^2 + 1/(2*h),0,0,0,...0;...
...
row n of the matrix:       0 0 0 0 0 0 ... r_n/h^2 - 1/(2*h) , -2*r_n/h^2, r_n/h^2 + 1/(2*h);...
row (n+1) of the matrix:   0 0 0 0 0 0 ... k/(2*h) , -alpha, -k/(2*h)]

Vector b: [Ta,0,0,0,0,....,0,-alpha*Tb]

If you want to stick to your arrangement, then do. I just saw at a quick glance that the rearranged equation for T_(n+1) cannot be correct because the terms T_(n+1) on the left-hand side and k* T_(n-1) on the right-hand side don't have the same unit.

[Casper Hae]

Ok I guess your method makes the equations "a little bit" easier.

But when plotting the solution this is what I get. And when changing N (number of intervals) the solution seem to change at the last boundary value (last vector in the y vector, if y = A\b). Is this true, because it seems quite weird to get pretty big differences in the solutions by varying N.

[Torsten]

Why is the plot not from 1 to 2 ?

Can you show your code ?

[Casper Hae]

Ta = 400; Tb = 22;

a = 1;

b = 2;

alpha = 10;

k = 65;

N = 50;

h = (b-a)/(N+1);

xp = (a+h:h:b-h)';

p = @(x) -2.*x/h^2;

q = @(x) x./h^2 + 1/(2*h);

r = @(x) x./h^2 - 1/(2*h);

1diag = p(xp);

2diag = q(xp);

3diag = r(xp);

A = diag(1diag) + diag(2diag(1:end-1),1) + diag(3diag(1:end-1),-1);

A(1,1) = 1;

A(1,2) = 0;

A(end,end-2) = k/(2*h);

A(end,end-1) = -alpha;

A(end,end) = -k/(2*h);

A = sparse(A);

bvec = zeros(N,1);

bvec(1) = Ta;

bvec(end) = -alpha*Tb;

y = A\bvec;

plot(x,y,"-o");

[Torsten]

Ta = 400; Tb = 22; 
a = 1; 
b = 2; 
alpha = 10;
k = 65;
N = 51;
h = (b-a)/(N-1);
xp = (a:h:b+h)';
np = numel(xp);
A = zeros(np,np);
A(1,1) = 1.0;
for i=2:np-1
  A(i,i-1) = xp(i)/h^2 - 1/(2*h);
  A(i,i) = -2*xp(i)/h^2;
  A(i,i+1) = xp(i)/h^2 + 1/(2*h);
end
A(end,end-2) = k/(2*h);
A(end,end-1) = -alpha;
A(end,end) = -k/(2*h);

B = zeros(np,1);
B(1) = Ta;
B(end) = -alpha*Tb;

sol = A\B;
plot(xp(1:end-1),sol(1:end-1))
hold on
c1 = 400;
c2 = -3780/65/(0.5+10/65*log(2));
yana = c1 + c2*log(xp);
plot(xp(1:end-1),yana(1:end-1))

[Torsten]

There was a slight mistake in the loop. I corrected it.

[Casper Hae]

Thanks for the help. But does this method has an order of accuracy of 2? I tried to evaluate the order of accuracy with your method and I believe it has an order of accuracy of 1. Or am I wrong?

[Torsten]

How did you check the order of accuracy "with my method" ?

[Casper Hae]

With this formula but substitute h with N (amount of intervals) as h and N are proportional to each other (h = (b-a)/(N-1)). p is the order of accuracy in this case.

[Torsten]

The quotient inside the log2(...) is a vector (where some components might be negative).

How do you take the log2 and result in a scalar value for this expression ?

[Casper Hae]

Creating a vector with the solutions for different N (N, N*2, N*4,N*8, etc.). Then evaluating with the formula above.

[Torsten]

Taking log2 of a vector gives a vector, not a scalar (p).

I don't want to lead you up the garden path, I just don't know how you come from a vector to a scalar (and the expression under the log2 might even have negative components so that log2 cannot be evaluated).

Do you mean perhaps

log2 (max(abs(u_h - u_h/2)) / max(abs(u_h/2-u_h/4)))

or something similar ?

[Casper Hae]

Creating a vector with solutions, e.g.

solutions = zeros(1,4);

if we have a function calculating T for different N (value(N)) we get

solutions(1) = value(10);

solutions(2) = value(20);

solutions(3) = value(40);

solutions(4) = value(80;

then applying the formula above we get

p = log2((solutions(1) - solutions(2))/(solutions(2)-solutions(3))

just as an example

[Torsten]

But then, you only evaluate the solution at one fixed r-coordianate. If you choose r=1 ,e.g, the error will always be 0 and you get order=Inf. Is this really legitimate ?

[Casper Hae]

Yeah let's say we check at r=1,5 with different N. How do you mean that the error will always be 0?

[Torsten]

At r=1, we have a fixed temperature that does not depend on h. So u_h = u_h/2 = u_h/4 = Ta. Thus p = log2(0/0) = NaN.

I think the definition of the order is more compliciated than you think it is. Maybe as I wrote: absolute value of maximum difference between the solutions for different h's.

[Casper Hae]

I've received p both with the formula above and by plotting it as a function of error. Both results give me an order of accuracy of 1. I got the equations from here. The absolute value of maxium difference is not really the order of accuracy from my understanding.

[Torsten]

|u_h - u | in this article is the local discretization error for the integration of an ODE.

It cannot be applied to the situation here.

[Casper Hae]

Okay. But by checking e.g.

solutions(4) - solutions(3)

and

solutions (3) - solutions(2)

and

solutions(2)-solutions(1)

we see that the convergence rate is linear. Doesn't this mean that the order of accuracy is 1 for this?

[Torsten]

We can continue discussion here, but the formula is not applicable for the situation here (at least not in the form you applied it).

The differential operators d^2T/dr^2 and dT/dr are discretized both with a spatial discretization scheme of order 2. Thus the total order of the method used is 2.

[Torsten]

In the theory of ordinary differential equations, the following relation between global error of a discretization and local discretization error holds:

global order of a discretization scheme = local discretization order - 1

Maybe this transfers to the situation here, although your comparisons are only based on differences in one point and thus cannot be seen as global error indicators.

I think

max(abs(u_h-u_h/2)) / max(abs(u_h/2-u_h/4))

is appropriate in this situation.

[Casper Hae]

Okay it makes sense. Where did you get that formula?

[Torsten]

http://eqworld.ipmnet.ru/en/methods/methods-ode/Niesen2004.pdf

Theorem 2.2

e.g. But it's the central result in the numerics of ordinary differential equations. It should be part of every book on this subject.

By the way: I could confirm your observation about the order. Also the expressions I used

max(abs(u_h-u_h/2)) / max(abs(u_h/2-u_h/4))

or

norm(u_h-u_h/2)/norm(u_h2-u_h/4)

are almost 2 in every case.

What I couldn't understand is that

max(u_h-u_analytic)/max(u_h/2-u_analytic) 

tends to 1.

I would have expected 2 also for this quotient.

[Torsten]

I found an error in my evaluation of the order of the method.

The terms

max(abs(u_h-u_h/2)) / max(abs(u_h/2-u_h/4))
norm(u_h-u_h/2)/norm(u_h2-u_h/4)
max(u_h-u_analytic)/max(u_h/2-u_analytic)

now all turned out to be near 4.

Thus the order of 2 of the method seems to be estabished.

I enclose the code - maybe you can search if there is still something wrong, but I believe no.

ntrials = 8; 
for j = 1:ntrials
    Ta = 400; 
    Tb = 22; 
    a = 1; 
    b = 2; 
    alpha = 10;
    k = 65;
    N = 2^(j-1)*50 + 1;  
    h = (b-a)/(N-1);
    xp = (a:h:b+h)';
    np = numel(xp);
    A = zeros(np,np);
    A(1,1) = 1.0;
    for i=2:np-1
        A(i,i-1) = xp(i)/h^2 - 1/(2*h);
        A(i,i) = -2*xp(i)/h^2;
        A(i,i+1) = xp(i)/h^2 + 1/(2*h);
    end
    A(end,end-2) = k/(2*h);
    A(end,end-1) = -alpha;
    A(end,end) = -k/(2*h);
    
    B = zeros(np,1);
    B(1) = Ta;
    B(end) = -alpha*Tb;
    
    A = sparse(A);
    B = sparse(B);
    
    sol = A\B;
    norm(A*sol-B)
    
    sol = sol(1:end-1);
    xp = xp(1:end-1);
    XP(:,j) = xp(1:2^(j-1):end);
    SOL(:,j) = sol(1:2^(j-1):end);   
end

diff1 = log2(max(abs((SOL(:,1:end-2)-SOL(:,2:end-1))))./max(abs((SOL(:,2:end-1)-SOL(:,3:end)))));
diff2 = log2(sqrt(sum((SOL(:,1:end-2)-SOL(:,2:end-1)).^2,1))./sqrt(sum((SOL(:,2:end-1)-SOL(:,3:end)).^2,1)));
c1 = Ta - log(a)*alpha*(Tb-Ta)/(k/b+alpha*log(b/a));
c2 = alpha*(Tb-Ta)/(k/b+alpha*log(b/a));
yana = c1 + c2*log(XP(:,1));
diff3 = log2(max(abs((SOL(:,1:end-1)-repmat(yana,1,size(SOL(:,1:end-1),2)))))./max(abs((SOL(:,2:end)-repmat(yana,1,size(SOL(:,2:end),2))))))
diff4 = log2(sqrt(sum((SOL(:,1:end-1)-repmat(yana,1,size(SOL(:,1:end-1),2))).^2,1))./sqrt(sum((SOL(:,2:end)-repmat(yana,1,size(SOL(:,2:end),2))).^2,1)))

figure(1)
plot(XP(:,1),[abs((SOL(:,1)-SOL(:,2))./(SOL(:,2)-SOL(:,3)))],"r",...
     XP(:,1),[abs((SOL(:,2)-SOL(:,3))./(SOL(:,3)-SOL(:,4)))],"g",...
     XP(:,1),[abs((SOL(:,3)-SOL(:,4))./(SOL(:,4)-SOL(:,5)))],"b",...
     XP(:,1),[abs((SOL(:,4)-SOL(:,5))./(SOL(:,5)-SOL(:,6)))],"m",...
     XP(:,1),[abs((SOL(:,5)-SOL(:,6))./(SOL(:,6)-SOL(:,7)))],"r",...
     XP(:,1),[abs((SOL(:,6)-SOL(:,7))./(SOL(:,7)-SOL(:,8)))],"g");
figure(2)    
plot(XP(:,1),[abs(SOL(:,1)-SOL(:,2))],"r",XP(:,1),[abs(SOL(:,2)-SOL(:,3))],"g",XP(:,1),[abs(SOL(:,3)-SOL(:,4))],"b",...
     XP(:,1),[abs(SOL(:,4)-SOL(:,5))],"r",XP(:,1),[abs(SOL(:,5)-SOL(:,6))],"g",XP(:,1),[abs(SOL(:,6)-SOL(:,7))],"b" ,...
     XP(:,1),[abs(SOL(:,7)-SOL(:,8))],"r")
figure(3)
plot(XP(:,1),[abs(SOL(:,1)-yana)],"r",XP(:,1),[abs(SOL(:,2)-yana)],"g",XP(:,1),[abs(SOL(:,3)-yana)],"b",XP(:,1),[abs(SOL(:,4)-yana)],"m",...
     XP(:,1),[abs(SOL(:,5)-yana)],"r",XP(:,1),[abs(SOL(:,6)-yana)],"g",XP(:,1),[abs(SOL(:,7)-yana)],"b",XP(:,1),[abs(SOL(:,8)-yana)],"m")
figure(4)
plot(XP(:,1),[SOL(:,1),yana],"r",XP(:,1),[SOL(:,2),yana],"g",XP(:,1),[SOL(:,3),yana],"b",XP(:,1),[SOL(:,4),yana],"m",...
     XP(:,1),[SOL(:,5),yana],"r",XP(:,1),[SOL(:,6),yana],"g",XP(:,1),[SOL(:,7),yana],"b",XP(:,1),[SOL(:,8),yana],"m")

Here is an example on how the order for finite-difference methods is usually computed (at least if the analytical solution is known). In the case it is not, your method - generalized to the treatment of a complete solution vector - seems adequate.

https://authors.library.caltech.edu/25061/4/NumMethChE84-Ch2-BVPforODEdiscretevar.pdf

page 75

[Torsten]

The overall accuracy of the method is 2 - also at the boundaries - since we discretized the boundary condition 2nd order and inserted it in a second-order discretization of the differential equation.

The evaluations of the results also show that the order of the method is 2 as you can readily see when you run the code from above. You must have made a mistake somewhere (like me at first).

If you want to work with the above equation to determine T_n, then use

xp = (a:h:b)';
np = numel(xp);
A = zeros(np,np);
A(1,1) = 1.0;
for i=2:np-1
    A(i,i-1) = xp(i)/h^2 - 1/(2*h);
    A(i,i) = -2*xp(i)/h^2;
    A(i,i+1) = xp(i)/h^2 + 1/(2*h);
end
A(end,end-2) = k/(2*h);
A(end,end-1) = -2*k/h;
A(end,end) = 3*k/(2*h) + alpha;

B = zeros(np,1);
B(1) = Ta;
B(end) = alpha*Tb;

This simply takes the equation

k*(T_(n-2) - 4*T_(n-1) + 3*T_n)/(2*h) = -alpha*(T_n - Tb)

as last equation to be solved in the linear system. Solution variable is T_n.

[Casper Hae]

Okay so I understand what you've done but have you used this discretization in the matrix:

k*(T_(n-2) - 4*T_(n-1) + 3*T_n)/(2*h) = -alpha*(T_n - Tb) ?

I don't really see where this discretization is used or how it is used since you approximated dT/dr with central differences instead of the discretization above (here). Or am I wrong? Thank you so much for this support btw.

[Torsten]

I don't really see where this discretization is used or how it is used since you approximated dT/dr with central differences instead of the discretization above (here). Or am I wrong?

This discretization in r=b

k*(T_(n-2) - 4*T_(n-1) + 3*T_n)/(2*h) = -alpha*(T_n - Tb)

is used in the piece of code of my last response:

A(end,end-2) = k/(2*h);
A(end,end-1) = -2*k/h;
A(end,end) = 3*k/(2*h) + alpha;

B(end) = alpha*Tb;

The complete code with this change would be

Ta = 400; 
Tb = 22; 
a = 1; 
b = 2; 
alpha = 10;
k = 65;
N = 51;
h = (b-a)/(N-1);
xp = (a:h:b)';
np = numel(xp);
A = zeros(np,np);
A(1,1) = 1.0;
for i=2:np-1
    A(i,i-1) = xp(i)/h^2 - 1/(2*h);
    A(i,i) = -2*xp(i)/h^2;
    A(i,i+1) = xp(i)/h^2 + 1/(2*h);
end
A(end,end-2) = k/(2*h);
A(end,end-1) = -2*k/h;
A(end,end) = 3*k/(2*h) + alpha;
B = zeros(np,1);
B(1) = Ta;
B(end) = alpha*Tb;
sol = A\B;
plot(xp,sol)
hold on
c1 = 400;
c2 = -3780/65/(0.5+10/65*log(2));
yana = c1 + c2*log(xp);
plot(xp,yana)

[Casper Hae]

Alright. So when discretizising with this discretization at r=b you don't have to account for eq (1) at all? Then it's so much easier lol.

Also, it is very weird that when I change to your last discretization in r=b the order of accuracy of the method suddenly drops to 1. It's like the new discretization has an order of accuracy of 1 (even though it says that it has 2)... do you know why?

[Torsten]

Also, it is very weird that when I change to your last discretization in r=b the order of accuracy of the method suddenly drops to 1. It's like the new discretization has an order of accuracy of 1 (even though it says that it has 2)... do you know why?

No, because the code i posted gives order 2 for both discretizations. Must be something wrong with your evaluation.

But you previously wrote that you got order 1 for the "old" discretization, too. Now you write that the order of accuracy "suddenly drops to 1". Did your result for the order of the "old" discretization change meanwhile ?