**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Crank-Nicholson method for solving telegraph quation

7 views (last 30 days)

Show older comments

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

##### 6 Comments

Torsten
on 30 Aug 2023

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

What exactly is your problem ?

HD
on 1 Sep 2023

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

### Answers (1)

Torsten
on 1 Sep 2023

Moved: Torsten
on 1 Sep 2023

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.

##### 25 Comments

HD
on 11 Sep 2023

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
on 11 Sep 2023

Edited: Torsten
on 11 Sep 2023

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.

Torsten
on 11 Sep 2023

Then have a read here:

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

Torsten
on 12 Sep 2023

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.

Torsten
on 12 Sep 2023

Edited: Torsten
on 12 Sep 2023

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
on 12 Sep 2023

Torsten
on 12 Sep 2023

Edited: Torsten
on 12 Sep 2023

HD
on 13 Sep 2023

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
on 13 Sep 2023

Edited: Torsten
on 13 Sep 2023

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:

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

HD
on 15 Sep 2023

Edited: Torsten
on 15 Sep 2023

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

HD
on 15 Sep 2023

Edited: HD
on 15 Sep 2023

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
on 15 Sep 2023

Torsten
on 15 Sep 2023

Edited: Torsten
on 15 Sep 2023

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

Torsten
on 16 Sep 2023

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
on 23 Sep 2023

Torsten
on 23 Sep 2023

Edited: Torsten
on 23 Sep 2023

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.

Torsten
on 23 Sep 2023

Edited: Torsten
on 23 Sep 2023

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

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)