Solve the partial differential equation using Crank-Nicolson method.
63 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hana Bachi
am 18 Feb. 2022
Kommentiert: Kuldeep Malik
am 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 Kommentare
Akzeptierte Antwort
Abolfazl Chaman Motlagh
am 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 Kommentare
Abolfazl Chaman Motlagh
am 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:
Kuldeep Malik
am 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
Weitere Antworten (1)
Siehe auch
Kategorien
Mehr zu General Applications finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!