Crank-Nicholson method for solving telegraph quation

Hello,

I'm trying to solve telegraph equation (transmission line) using Crank-Nicholson in MATLAB, but I'm stuck with Crank-Nicholson difference scheme. Can someone help?

I'm using simplified model

[Torsten]

The differencing scheme can be found everywhere in the internet, e.g. here:

https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method

What exactly is your problem ?

[HD]

I'dont know how to apply scheme to first element in equation

[Torsten]

By writing your equation as a system of two first-order equations

du/dt = v
dv/dt = d^2u/dx^2 / LC

or

du/dt + 1/sqrt(LC)*du/dx = v
dv/dt - 1/sqrt(LC)*dv/dx = 0 

https://en.wikipedia.org/wiki/Wave_equation

[HD]

Is this correct

[Torsten]

You can try this, but it's not Crank-Nicolson. I'd call it Explicit Euler.

[HD]

So can i write it like this

My question now is can I substitute in (2) using (1) and how can I write ?

If i use this eq , and if i write what is correct way to write ? Thanks

Answer 1

[Torsten]

It's unusual that the lower index denotes the time discretization - thus I will write y_{i}^{n} for y at time t(n) in grid point x(i).

Your CN discretization reads

(u_{i}^{n+1}-u_{i}^{n})/dt = (v_{i}^(n+1)+v_{i}^{n})/2
(v_{i}^{n+1}-v_{i}^{n})/dt = 1/(LC) * (u_{i+1}^{n+1}-2*u_{i}^{n+1}+u_{i-1}^{n+1} + u_{i+1}^{n}-2*u_{i}^{n}+u_{i-1}^{n})/(2*dx^2)

or

u_{i}^{n+1}/dt - v_{i}^(n+1)/2 = u_{i}^{n}/dt +v_{i}^{n}/2
-1/(LC)*(u_{i-1}^{n+1}-2*u_{i}^{n+1}+u_{i+1}^{n+1})/(2*dx^2) + v_{i}^{n+1}/dt = v_{i}^{n}/dt + 1/LC*(u_{i+1}^{n}-2*u_{i}^{n}+u_{i-1}^{n})/(2*dx^2)

This is a system of linear equations in the unknowns u_{i}^{n+1} and v_{i}^{n+1} (2<=i<=nx-1).

For indices i = 1 and i = nx, you have to incorporate the spatial boundary conditions for u and v = du/dt.

[HD]

Hello again, this is my code. I wrote and put that in second equation. Is this correct way?

%% Define parameters

a = 0;

b = 1;

L = 1.0; %Length of the wire (1 meter)

nx = 10; %Number of spatial grid points (including boundaries)

nt = 30; %Number of time steps

dx = (b-a) / (nx - 1); % Spatial step

t0 = 0;

tf = 2 ; % Final time

dt = (tf - t0)/(nt-1); % Time step

LC = 1 ; % L/C constant (adjust as needed)

x = a:dx:b;

t = t0:dt:tf;

% Time-stepping loop

u = zeros(nx, nt);

v = zeros(nx, nt);

u(:,1) = sin(pi.*x); % Sinusoidal voltage source

v(:,1) = pi*sin(pi.*x); % Time derivative of voltage

u(:,2) = sin(pi*x)*(1+dt);

v(:,2) = pi*sin(pi.*x)*(1+dt);

for n = 2:nt-1

for i = 2:nx-1

u(i, n+1) =((dt/2*dx^2*LC)*(u(i+1,n+1)+u(i-1,n+1)+u(i+1,n)-2*u(i,n)+u(i-1,n))+2*u(i,n)/dt+2*v(i,n))/((2/dt)+(dt/(dx^2*LC)));

v(i, n+1) = (2/dt)*(u(i,n+1)-u(i,n))-v(i,n);

end

end

[Torsten]

Is this correct way?

No. The terms on the left-hand side are the unknowns, the terms on the right-hand side are known. So

u_{i}^{n+1}/dt - v_{i}^(n+1)/2 = u_{i}^{n}/dt +v_{i}^{n}/2
-1/(LC)*(u_{i-1}^{n+1}-2*u_{i}^{n+1}+u_{i+1}^{n+1})/(2*dx^2) + v_{i}^{n+1}/dt = v_{i}^{n}/dt + 1/LC*(u_{i+1}^{n}-2*u_{i}^{n}+u_{i-1}^{n})/(2*dx^2)

can be written as

A*[u^(n+1);v^(n+1)] = b(u^(n),v^(n))

Now solve for [u^(n+1);v^(n+1)] as

[u^(n+1);v^(n+1)] = A\b(u^(n),v^(n))

And take care of the boundary values for u and v. You didn't treat them in your code because your loops run from 2 to nx-1.

[HD]

I can't figure out how you got this, sorry

[Torsten]

Then have a read here:

https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method

The two 1d-examples should show you how to derive A and b for your case.

[HD]

Is this now correct for i = 1,2, 3?

[Torsten]

I suggest you order your variables as (u_1,v_1,u_2,v_2,...) instead of (u_1,u_2,...,v_1,v_2,...) because this will give a small bandwidth for A and thus a more efficient solving of A*x = b.

So, if you order your variables as [u_1,v_1,u_2,v_2,u_3,v_3,u_4,v_4,u_5,v_5] and u_1, v_1, u_5, v_5 are the boundary values, the matrix A looks like

A = [1/dt -1/2 0 0 0 0 0 0 0 0;
     * * * * * * * * * *
     0 0 1/dt -1/2 0 0 0 0 0 0;
     -r 0 2*r 1/dt -r 0 0 0 0 0;
     0 0 0 0 1/dt -1/2 0 0 0 0;
     0 0 -r 0 2*r 1/dt -r 0 0 0;
     0 0 0 0 0 0 1/dt -1/2 0 0;
     0 0 0 0 -r 0 2*r 1/dt -r 0;
     0 0 0 0 0 0 0 0 1/dt -1/2;
     * * * * * * * * * *]    

with

r = 1/(LC*2*dx^2)

.

The rows indicated with "*" must be filled with the boundary conditions.

[HD]

My system looks like this

Is this ok?

[Torsten]

Is this ok?

No. The matrix A in your notation must read

[1/dt 0 0 0 -1/2 0 0 0;
 * * * * * * * *;
 0 1/dt 0 0 0 -1/2 0 0;
 -r 2r -r 0 0 1/dt 0 0;
 0 0 1/dt 0 0 0 -1/2 0;
 0 -r 2r -r 0 0 0 1/dt;
 0 0 0 1/dt 0 0 0 -1/2;
 * * * * * * * *]

with

r = 1/(LC*2*dx^2)

and the vector you solve for is

[u_0^n+1,u_1^(n+1),u_2^(n+1),u_3^(n+1),v_0^(n+1),v_1^(n+1),v_2^(n+1),v_3^(n+1)]

The eight stars in rows 2 and 8 stand for the coefficients of the boundary condition for u_0^(n+1) and u_3^(n+1).

I cannot understand why it's so difficult for you to put the two equations for each grid point I gave you into a matrix form. After making an attempt, just multiply out with pencil and paper to see if you reproduce the equations given.

[HD]

Thank you for your time and patience, I will try the way you explain to me. It is difficult to me because i don’t understand boundary conditions and how to represent system like that in matrix form.

[Torsten]

I also don't understand why you program the numerical method (Crank-Nicholson) on your own if you have so little experience with numerical methods. Using ODE15s, you only need to supply the time derivatives for u and v in the grid points, and the solver does everything else for you. Is it a challenge or a homework assignment you are working on ?

[HD]

This is homework assignement for me and it's very important. Sorry for my experience and I appreciate your effort and time. I'm last year student and this method is assigned to me to solve telegraph equation (transmission line).

[HD]

Can you please help me to solve this, I don't know how to put boundary conditions. Equation is with , , , . If my system is now correct?

Thanks in advance.

[Torsten]

Looks fine to me. But why is the c*du/dt - term missing in your original equation ? As written, you solve the wave equation, not the telegraph equation:

https://en.wikipedia.org/wiki/Telegrapher%27s_equations

Concerning your questions:

Initialize

[u_1^0,v_1^0,u_2^0,v_2^0,u_3^0,v_3^0,u_4^0,v_4^0,u_5^0,v_5^0] = [0,0,sin(pi/4),sin(pi/4),sin(pi/2),sin(pi/2),sin(3/4*pi),sin(3/4*pi),sin(pi),sin(pi)]

choose the second row on both sides as

[1 0 0 0 0 0 0 0 0 0],

(meaning u(0) = sin(pi) for all times) and the tenth row on both sides as

[0 0 0 0 0 0 0 0 1 0]

(meaning u(1) = sin(pi) for all times) and start solving for

[u_1^1,v_1^1,u_2^1,v_2^1,u_3^1,v_3^1,u_4^1,v_4^1,u_5^1,v_5^1]

and so on.

Note that you need to factorize the matrix on the left-hand side only once and then solve the system for different right-hand sides. Read the chapter

Solving for Several Right-Hand Sides

under

https://uk.mathworks.com/help/matlab/math/systems-of-linear-equations.html

[HD]

Hello again, hope I get it right

% Parameters to define the heat equation and the range in space and time
L = 1.; % Lenth of the wire
T =1.; % Final time
maxk = 10; % Number of time steps
dt = T/maxk;
n = 4.; % Number of space steps
dx = L/n;
LC = 1;
r = 1/(1*2*LC*dx^2);
%%
A=[ 1/dt -0.5 0 0 0 0 0 0 0 0;
    1 0 0 0 0 0 0 0 0 0;
    0 0 1/dt -0.5 0 0 0 0 0 0;
   -r 0 2*r 1/dt -r 0 0 0 0 0;
    0 0 0 0 1/dt -0.5 0 0 0 0;
    0 0 -r 0 2*r 1/dt -r 0 0 0;
    0 0 0 0 0 0 1/dt -0.5 0 0;
    0 0 0 0 -r 0 2*r 1/dt -r 0 ;
    0 0 0 0 0 0 0 0 1/dt -0.5;
    0 0 0 0 0 0 0 0 1 0];
dA = decomposition(A);
%% 
B = [1/dt 0.5 0 0 0 0 0 0 0 0;
    1 0 0 0 0 0 0 0 0 0;
    0 0 1/dt 0.5 0 0 0 0 0 0;
    r 0 -2*r 1/dt r 0 0 0 0 0;
    0 0 0 0 1/dt 0.5 0 0 0 0;
    0 0 r 0 -2*r 1/dt r 0 0 0;
    0 0 0 0 0 0 1/dt 0.5 0 0;
    0 0 0 0 r 0 -2*r 1/dt r 0 ;
    0 0 0 0 0 0 0 0 1/dt 0.5;
    0 0 0 0 0 0 0 0 1 0];
%% Boundary conditions
values = [0, 0, sin(pi/4), sin(pi/4), sin(pi/2), sin(pi/2), sin(3/4*pi),sin(3/4*pi), sin(pi),sin(pi)];
u(:,1) = values';
%%
for i = 1:maxk
  u(:,i+1) = dA\(B*u(:,i));
end; u

u = 10×11

0 0 0 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 0 0 0 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.7071 0.7438 0.7124 0.6158 0.4627 0.2673 0.0474 -0.1768 -0.3849 -0.5577 -0.6794 0.7071 0.0272 -0.6553 -1.2777 -1.7831 -2.1252 -2.2727 -2.2121 -1.9488 -1.5071 -0.9274 1.0000 1.0519 1.0075 0.8708 0.6544 0.3780 0.0671 -0.2501 -0.5443 -0.7887 -0.9608 1.0000 0.0384 -0.9267 -1.8069 -2.5217 -3.0055 -3.2141 -3.1283 -2.7561 -2.1314 -1.3116 0.7071 0.7438 0.7124 0.6158 0.4627 0.2673 0.0474 -0.1768 -0.3849 -0.5577 -0.6794 0.7071 0.0272 -0.6553 -1.2777 -1.7831 -2.1252 -2.2727 -2.2121 -1.9488 -1.5071 -0.9274 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000

plot(0:dx:L,u(1:2:2*(n+1),:))

[Torsten]

To get the solution at the final time, the loop must run up to 10 instead of 9:

for i = 1:maxk
    u(:,i+1) = dA\(B*u(:,i));
end

Do you have a solution to compare with ?

At least you should finally plot u :

plot(0:dx:L,u(1:2:2*(n+1),:))

[HD]

I have analytic solution

a = 0;
b = 1;
L = 1.0;             %Length of the wire (1 meter)
nx = 5;           %Number of spatial grid points (including boundaries)
nt = 10;          %Number of time steps
dx = (b-a) / (nx - 1);  % Spatial step
t0 = 0;
tf = 1 ;  % Final time
dt = (tf - t0)/(nt);        % Time step

LC = 1.0 ;          % L/C constant (adjust as needed)

x = a:dx:b;
t = t0:dt:tf;
s=dt^2/dx^2;
UA = zeros(nx,nt);
for i = 1:nx
	for j = 1:nt
	UA(i,j) = sin(pi*x(i))*(cos(pi*t(j))+sin(pi*t(j))/pi);
    end
end
plot(0:dx:L,UA)

[Torsten]

And do the curves coincide for both analytical and numerical solution if you plot them together ? (I think nx must be 5 in the code above).

[HD]

I changed nx to 5 and plot both solutions. Dashed curves represent analytical solution.

[Torsten]

And what is your conclusion ?

I'd say a code that automatically generates the matrices on both sides of your linear system depending on the number of grid points would be helpful. More grid points - better precision. But it looks promising.

But I'm surprised that the numerical solution looks different from the one I posted above (which came from your code).

[HD]

What is different? I will increase grid points, this is just for test.

[Torsten]

What is different?

I think the number of plotted curves is different.

% Parameters to define the heat equation and the range in space and time
L = 1.; % Lenth of the wire
T =1.; % Final time
maxk = 10; % Number of time steps
dt = T/maxk;
n = 4.; % Number of space steps
dx = L/n;
LC = 1;
r = 1/(1*2*LC*dx^2);
%%
A=[ 1/dt -0.5 0 0 0 0 0 0 0 0;
    1 0 0 0 0 0 0 0 0 0;
    0 0 1/dt -0.5 0 0 0 0 0 0;
   -r 0 2*r 1/dt -r 0 0 0 0 0;
    0 0 0 0 1/dt -0.5 0 0 0 0;
    0 0 -r 0 2*r 1/dt -r 0 0 0;
    0 0 0 0 0 0 1/dt -0.5 0 0;
    0 0 0 0 -r 0 2*r 1/dt -r 0 ;
    0 0 0 0 0 0 0 0 1/dt -0.5;
    0 0 0 0 0 0 0 0 1 0];
dA = decomposition(A);
%% 
B = [1/dt 0.5 0 0 0 0 0 0 0 0;
    1 0 0 0 0 0 0 0 0 0;
    0 0 1/dt 0.5 0 0 0 0 0 0;
    r 0 -2*r 1/dt r 0 0 0 0 0;
    0 0 0 0 1/dt 0.5 0 0 0 0;
    0 0 r 0 -2*r 1/dt r 0 0 0;
    0 0 0 0 0 0 1/dt 0.5 0 0;
    0 0 0 0 r 0 -2*r 1/dt r 0 ;
    0 0 0 0 0 0 0 0 1/dt 0.5;
    0 0 0 0 0 0 0 0 1 0];
%% Boundary conditions
values = [0, 0, sin(pi/4), sin(pi/4), sin(pi/2), sin(pi/2), sin(3/4*pi),sin(3/4*pi), sin(pi),sin(pi)];
u(:,1) = values';
%%
for i = 1:maxk
  u(:,i+1) = dA\(B*u(:,i));
end

x = 0:dx:L;
t = 0:dt:T;
for i = 1:numel(x)
  for j = 1:numel(t)
      ua(i,j) = sin(pi*x(i))*(cos(pi*t(j))+sin(pi*t(j))/pi);
  end
end

[r,c] = size(ua);
cc = lines(c);

hold on
for j = 1:c
  plot(x,u(1:2:2*(n+1),j),'-','Color',cc(j,:))
  plot(x,ua(1:n+1,j),'o','Color',cc(j,:))
end
hold off

[HD]

Hello this is great solution, thank you very much. Can you please just help me with one more example where boundary condition

[Torsten]

Say u(1,t) = f(t). Then v(1,t) = f'(t).

Thus Crank-Nicolson says that you should implement it as

1/2*(u_5^(n+1) + u_5^n) = 1/2*( f(t^(n+1)) + f(t^n))
1/2*(v_5^(n+1) + v_5^n) = 1/2*( f'(t^(n+1)) + f'(t^n))

or

1/2*u_5^(n+1) = -1/2*u_5^n + 1/2*( f(t^(n+1)) + f(t^n))
1/2*v_5^(n+1) = -1/2*v_5^n + 1/2*( f'(t^(n+1)) + f'(t^n))

This shows you how you have to modify the last (2x2) matrix in A and B and that you have to add a vector b on the right-hand side of your iteration given by

b = [0; 0; 0; 0; 0; ... ;1/2*( f(t^(n+1)) + f(t^n));1/2*( f'(t^(n+1)) + f'(t^n))]   

By the way: You can of course use this method also if u(1,t) = 0 (as in the reference case upto now).

Be careful if u(1,0) or v(1,0) are not equal to your initial condition for u and v at x = 1.

[HD]

Thank you very much

[Torsten]

I'm not sure whether this averaging 1/2*(unew + uold) in Crank-Nicolson also applies to algebraic equations like the boundary conditions.

You might want to compare the results to simply setting

u_5^(n+1) = f(t^(n+1)) 
v_5^(n+1) = f'(t^(n+1)) 

thus setting the (2x2) matrix in A to [1 0;0 1], in B to [0 0 ; 0 0] and the vector b to

b = [0; 0; 0; 0; 0; ... ; f(t^(n+1)) ;f'(t^(n+1)) ]