The differencing scheme can be found everywhere in the internet, e.g. here:
What exactly is your problem ?
Crank-Nicholson method for solving telegraph quation
12 次查看（过去 30 天）
显示 更早的评论
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 个评论
HD
2023-9-1
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 回答（1 个）
Torsten
2023-9-1
移动：Torsten
2023-9-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 个评论
HD
2023-9-11
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
Torsten
2023-9-11
编辑：Torsten
2023-9-11
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
2023-9-11
Then have a read here:
The two 1d-examples should show you how to derive A and b for your case.
Torsten
2023-9-12
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
2023-9-12
编辑：Torsten
2023-9-12
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
2023-9-12
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
2023-9-12
编辑：Torsten
2023-9-12
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
2023-9-13
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.
Torsten
2023-9-13
编辑：Torsten
2023-9-13
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
2023-9-15
编辑：Torsten
2023-9-15
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
2023-9-15
编辑：HD
2023-9-15
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
2023-9-15
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
2023-9-15
编辑：Torsten
2023-9-15
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
2023-9-16
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
2023-9-23
Hello this is great solution, thank you very much. Can you please just help me with one more example where boundary condition 
Torsten
2023-9-23
编辑：Torsten
2023-9-23
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
2023-9-23
编辑：Torsten
2023-9-23
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)) ]
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)