Crank-Nicholson method for solving telegraph quation
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 31 Aug 2023
I'dont know how to apply scheme to first element in equation
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
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 
?
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
, and if i write what is correct way to write ? Thanks Answers (1)
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?
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
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
on 11 Sep 2023
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.
HD
on 12 Sep 2023
My system looks like this 
Is this ok?
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
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
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?
with , , , . If my system is now correct?
Thanks in advance.
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
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
plot(0:dx:L,u(1:2:2*(n+1),:))
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),:))
Torsten
on 15 Sep 2023
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).
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
on 15 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
Hello this is great solution, thank you very much. Can you please just help me with one more example where boundary condition 
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
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)) ]
Categories
Find more on Calendar in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
