Solve the partial differential equation using Crank-Nicolson method.
50 views (last 30 days)
Show older comments
Hana Bachi
on 18 Feb 2022
Commented: Kuldeep Malik
on 10 Aug 2023
Hello everyone,
you may see my questions here from time to time. I'm a beginner in Mathlab and trying to solve some equations, then ask for advice if I accounter any errors.
I tried to solve a PDE using Crank-Nicolson method. I wrote a code but I'm still getting some errors, which didn't allow me to get the final results. any advice, please.
I attached the code and the question.
Thank you
0 Comments
Accepted Answer
Abolfazl Chaman Motlagh
on 18 Feb 2022
Hi.
you forget to put index in x vector in line 70, in equation you calculate Ur. that's the error.
Ur(i,j) = 0.5*erfc(0.5*(x(i)-t(j+1))./(B*sqrt(t(j+1))))+0.5*exp(x(i)/B^2).*erfc(0.5*(x(i)+t(j+1))/(B*sqrt(t(j+1))));
and also you put j+1 in t, but j is from 1 to M=1651. but at j=1651 the t(1652) doesn't exist. so correct that too.
8 Comments
Abolfazl Chaman Motlagh
on 21 Feb 2022
Hi hana.
i check your code, i couldn't figure it out exactly. there were a lot of ambiguities for me. specially variable d that you change it in every loop but never use it!
i tried to solve it myself. since my first try failed, the problem really got me :) . i literally solve the problem (discretization) more than 5 times from scratch. finally it works. the linear system is really sensitive to one little mistake in calculations.
i attached my code and also my formula for solution. hopes it is what you need.
here's final solution for 2 different discretization and the analytical approximation the pdf provide:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/902230/image.jpeg)
Kuldeep Malik
on 10 Aug 2023
Hello Dear
Can you help me with the code
I know something is wrong with the right boundary condition.
I have used mittag leffler function for initial and boundry conditions.
% clear; clc;
dx = 1/1000;
dt = 1/300000; % 1/3300
i_max = round(1/dx) + 1;
n_min = 1;
n_max = 10;
a=0.5;
La=(dx/dt)^a;
%beta = sqrt(1/878);
%r = dt / (4*dx);
%alpha = (dt*(beta^2))/(2*(dx^2));
% a = r-alpha;
b = 2*La;
%c = - (r+alpha);
%b_ = 1 - 2*alpha;
A = b * eye(i_max-2);
for k=1:size(A,1)
if k<size(A,1)
A(k,k+1) = 1;
end
if k>1
A(k,k-1) = -1;
end
end
% A(end) = A(end) + a;% right boundry condition
% u(nx,j)=((mlf(a,1,((nx-1).*dx).^a,7)-mlf(a,1,-1.*(((nx-1).*dx).^a),7))./2).*((mlf(a,1,((j-1).*dt).^a,7)+mlf(a,1,-1.*((j-1).*dt).^a,7))./2)-((mlf(a,1,((nx-1).*dx).^a,7)+mlf(a,1,-1.*(((nx-1).*dx).^a),7))./2).*((mlf(a,1,((j-1).*dt).^a,7)-mlf(a,1,-1.*((j-1).*dt).^a,7))./2);
A;
U = zeros(i_max,n_max);
% U(:,1) = 0;
% initial Condition
for i=2:i_max
U(i,1)=(mlf(a,1,((i-1).*dx).^a,7)-mlf(a,1,-1.*(((i-1).*dx).^a),7))./2;
end
% U(1,:) = 1;
% left and right boundry condition
for j=1:n_max
U(1,j)=-((mlf(a,1,((j-1).*dt).^a,7)-mlf(a,1,-1.*(((j-1).*dt).^a),7))./2); %left
U(i_max,j)=((mlf(a,1,((i_max-1).*dx).^a,7)-mlf(a,1,-1.*(((i_max-1).*dx).^a),7))./2).*((mlf(a,1,((j-1).*dt).^a,7)...
+mlf(a,1,-1.*((j-1).*dt).^a,7))./2)-((mlf(a,1,((i_max-1).*dx).^a,7)+mlf(a,1,-1.*(((i_max-1).*dx).^a),7))./2).*((mlf(a,1,((j-1).*dt).^a,7)...
-mlf(a,1,-1.*((j-1).*dt).^a,7))./2); %right
end
U;
for n = n_min:n_max-1
RHS = zeros(i_max-2,1);
for l = 1:size(A,1)
if l==1
RHS(l) = b * U(l+1,n) - U(l+2,n) + U(l,n) + U(1,n+1);
elseif l==size(A,1)
RHS(l) = b * U(l+1,n)- U(l+2,n) + U(l,n) - U(l+2,n+1);
else
RHS(l)=b * U(l+1,n) - U(l+2,n) + U(1,n);
end
end
u_local = A\RHS;
U(2:end-1,n+1) = u_local;
U(end,n+1) = u_local(end); % right boundry condition
end
U;
Ur = zeros(i_max,n_max);
Error = zeros(i_max,n_max);
x = @(i) ((i-1)*dx);
t = @(i) ((i-1)*dt);
%Exact Solution is here and Error
for jjj=1:n_max
for iii=1:i_max
Ur(iii,jjj)=((mlf(a,1,((iii-1).*dx).^a,7)-mlf(a,1,-1.*(((iii-1).*dx).^a),7))./2).*((mlf(a,1,((jjj-1).*dt).^a,7)+mlf(a,1,-1.*((jjj-1).*dt).^a,7))./2)-((mlf(a,1,((iii-1).*dx).^a,7)+mlf(a,1,-1.*(((iii-1).*dx).^a),7))./2).*((mlf(a,1,((jjj-1).*dt).^a,7)-mlf(a,1,-1.*((jjj-1).*dt).^a,7))./2);
Error(iii,jjj)=abs(Ur(iii,jjj)-U(iii,jjj));
end
end
U
Ur
Error
% for j=1:n_max %Time Loop
% for i=1:i_max %Space Loop
% Ur(i,j) = 0.5*erfc(0.5*(x(i)-t(j))/(beta*sqrt(t(j))))+0.5*exp(x(i)/(beta^2)).*erfc(0.5*(x(i)+t(j))/(beta*sqrt(t(j))));
% Error(i,j)=abs(Ur(i,j)-U(i,j));
% end
% end
% %% Final Plot
% x = 0:dx:1;
% figure;
% plot(x,U(:,end),'LineWidth',2); hold on;grid on;
% plot(x,Ur(:,end),'LineWidth',2);
% title(['T = ' num2str(dt * n_max)]);
% pause(0.5);
% %% Animation over Time
% x = 0:dx:1;
% figure;
% for i=1:size(U,2)
% plot(x,U(:,i),'LineWidth',2); hold on;grid on;
% plot(x,Ur(:,i),'LineWidth',2);
% title(num2str((i-1)*dt));
% hold off;
% pause(dt)
% end
More Answers (1)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!