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
9 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 ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1468891/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1468891/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1470766/image.png)
My question now is can I substitute
in (2) using (1) and how can I write
?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1470771/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1470776/image.png)
If i use this eq
, and if i write
what is correct way to write
? Thanks
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1470781/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1470786/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1470791/image.png)
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?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1479251/image.png)
%% 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
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
on 12 Sep 2023
Edited: Torsten
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?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1481756/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1481761/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1481766/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1481771/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1481776/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1481781/image.png)
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),:))
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1483701/image.png)
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
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).
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1484247/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1490947/image.png)
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
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 (한국어)