Solve the partial differential equation using Crank-Nicolson method.

73 次查看（过去 30 天）
Hana Bachi
Hana Bachi 2022-2-18
评论： Kuldeep Malik 2023-8-10
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

请先登录，再回答此问题。

采纳的回答

Abolfazl Chaman Motlagh
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 个评论
显示 6更早的评论
Abolfazl Chaman Motlagh
Abolfazl Chaman Motlagh 2022-2-21
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
Kuldeep Malik 2023-8-10
@Abolfazl Chaman Motlagh
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

请先登录，再进行评论。

更多回答（1 个）

Hana Bachi
Hana Bachi 2022-2-22
Hello Abolfazi,
thank you so much for giving from your time to solve this problem, I deeply appreciate it.
The d term in my code is RHS in yours.
I found that I forget to multiply some terms by B^2, and I added dt/2 somewhere to, this is why it was showing big difference between the two solution.
But I think yours make more sence then the one I got
here is mine after writing the code

类别

在 Help Center 和 File Exchange 中查找有关 General Applications 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by