- #1
Here is the codeDrClaude said:There is a bug on line 42
If you don't post the code, there is not much anyone can do.
%% Crank-Nicolson Method
clear variables
close all
% 1. Space steps
xa = 0;
xb = 1;
dx = 1/40;
N = (xb-xa)/dx ;
x = xa:dx:xb;
%2.Time steps
ta = 0;
tb = 0.5;
dt = 1/3300;
M = (tb-ta)/dt ;
t = ta:dt:tb;
%3. Controling Parameters
%4. Define equations A, B , C and phi(x,t)
A = @(x,t) (50/3)*(x-0.5+4.95*t);
B = @(x,t) (250/3)*(x-0.5+0.75*t);
C = @(x,t) (500/3)*(x-0.375);
phi = @(x,t)(0.1*exp(-A(x,t)) +0.5*exp(-B(x,t))+exp(-C(x,t)))./...
(exp(-A(x,t)) +exp(-B(x,t))+exp(-C(x,t)));
% 5 Initial and boundary conditions
f = @(x) phi(x,t(1)); % Initial condition
g1 = @(t)phi(x(1),t); % Left boundary condition
g2 = @(t)phi(x(N+1),t); % Right boundary condition
r = dt / (2*dx^2);
a = -0.003*r;
b = 1 + 2*0.003*r;
c = -r*0.003;
d =(1+0.003*2*r);
e =1+0.003*2*r;
h =0.003*r;
% 6 Implementation of the explicit method
U = zeros(N+1,M+1);
U(2:N,1) = f(x(2:N)); % Put in the initial condition
U(1,:) = g1(t); % The boundary conditions, g1 and g2 at x = 0 and x = 1
U(N+1,:) = g2(t);
for i = 2:N
RHS= (1+2*0.003*r)*U(i,N)+0.003*r*U(i+1,N)+0.003*r*U(i-1,N)-r*dx*U(i,N)*(U(i,N)-U(i-1,N));
end
for j=2:M % Time Loop
for i= 2:N % Space Loop
RHS= (1+2*0.003*r)*U(i,N+1)+0.003*r*U(i+1,N+1)+0.003*r*U(i-1,N+1)-r*dx*U(i,N+1)*(U(i,N+1)-U(i-1,N+1));
end
end
% Make some plots
T= 0:.1:M;
V = [];
for i= 1: length(T)
P = find(t==T(i));
V = [V P];
end
figure
subplot(131)
for j = 1:length(V)
hold on
plot(x,U(:,V(j)),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))
end
legend('-DynamicLegend','location','bestoutside');
a = ylabel('U-Values');
set(a,'Fontsize',14);
a = xlabel('X-Values');
set(a,'Fontsize',14);
a=title('Crank-Nicolson Solution');
set(a,'Fontsize',16);
grid;
% disp(u(:,V)'); % Each row corresponds to a particular value of t and
% Each column corresponds to a particular value of x
% Implement the exact solution and compare it to the exact solution
Exact =@(x,t) phi(x,t);
subplot(132)
for j = 1:length(V)
hold on
plot(x,Exact(x,t(V(j))),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))
end
legend('-DynamicLegend','location','bestoutside');
a = ylabel('U-Values');
set(a,'Fontsize',14);
a = xlabel('X-Values');
set(a,'Fontsize',14);
a=title(' Analytical Solution');
set(a,'Fontsize',16);
grid;
[X,T] = meshgrid(x,t);
% Exact2 = (0.1*exp(-A(X,T)) +0.5*exp(-B(X,T))+exp(-C(X,T)))./...
% (exp(-A(X,T)) +exp(-B(X,T))+exp(-C(X,T)));
Error =abs(Exact(X,T)'-U);
subplot(133)
for j = 1:length(V)
hold on
plot(x,Error(:,(V(j))),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))
end
legend('-DynamicLegend','location','bestoutside');
a = ylabel('Error');
set(a,'Fontsize',14);
a = xlabel('X-Values');
set(a,'Fontsize',14);
a=title(' Absolute Error ');
set(a,'Fontsize',16);
grid;
I used this dt = 1/3300pasmith said:What's the maximum timestep you can use and still have the method be stable? Are you exceeding it?
This is linear. I assume 0U is a typo though?hanabachi said:actually the equation is
dU/dt=0U*du/dx+0.003*dU2/dx2
Actually the partial differential equation is:hanabachi said:actually the equation is
dU/dt=0U*du/dx+0.003*dU2/dx2
0 should be negative signs (-), so the equation is the multiplication of U by its first derivative dU/dx. So the equation is not linearjoshmccraney said:This is linear. I assume 0U is a typo though?
bigfooted said:Actually the partial differential equation is:
$$ \frac{\partial u}{\partial t} = -u\frac{\partial u}{\partial x} + 0.003 \frac{\partial^2 u }{\partial x^2}, x \in (0,1), t>0$$
But the question is to solve it using the Crank-Nicolson method. Do you have the difference scheme for this PDE and can you reproduce it here? We can continue from there.
Please be aware that a lot of people simply refuse to download data trash, i.e. documents that have to be deleted after use. I, for example. If it is too troublesome for you to type it out here, then it is too troublesome for me to download whatever you want to force me to.hanabachi said:Yes, this is what I'm using to solve the pde. please check the attached file
yes, here it isbigfooted said:Can you please type it here in latex? My eyes are still hurting from looking at this document :-)
bigfooted said:Actually the partial differential equation is:
$$ \frac{\partial u}{\partial t} = -u\frac{\partial u}{\partial x} + 0.003 \frac{\partial^2 u }{\partial x^2}, x \in (0,1), t>0$$
But the question is to solve it using the Crank-Nicolson method. Do you have the difference scheme for this PDE and can you reproduce it here? We can continue from there.
Yes, that what it should looks likebigfooted said:OK, nice. So the left-hand side creates the matrix A and the right hand side created the vector b. Can you show now how to construct a row in matrix A and vector b? So what does the first and second row look like for instance?
row 1: $$[A_{11} A_{21} ... A_{N1}] \cdot U_1 = b_1$$
The Crank Nicolson method is a numerical method used to solve partial differential equations (PDEs). It is a combination of the explicit Euler method and the implicit Euler method, and is known for its stability and accuracy.
The Crank Nicolson method uses a finite difference approximation to discretize the PDE and then solves the resulting system of equations using a tridiagonal matrix algorithm. This method takes into account both the current and future time steps, resulting in a more accurate solution.
The Crank Nicolson method can be used to solve a wide range of PDEs, including parabolic, hyperbolic, and elliptic equations. It is particularly useful for solving diffusion and heat transfer problems.
No, the accuracy of the Crank Nicolson method depends on the specific PDE being solved and the chosen time and space discretization. In some cases, other methods may be more accurate. It is important to carefully consider the problem at hand and choose the most appropriate method.
The Crank Nicolson method has several advantages, including its stability, accuracy, and ability to handle a wide range of PDEs. It is also relatively easy to implement and can handle non-uniform grids. Additionally, it is a second-order method, meaning that it converges to the exact solution at a faster rate than first-order methods.