Solve the partial differential equation using Crank-Nicolson method.

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

Answer 1

[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.

[Hana Bachi]

Hi Abolfazi,

I'm still getting this :

Index exceeds the number of array elements. Index must not exceed 1651.

Error in Cp32nd (line 70)

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))));

[Abolfazl Chaman Motlagh]

yeah exactly, i warn you about that in last line. you can do this:

Ur(i,j) = 0.5*erfc(0.5*(x(i)-t(j))./(B*sqrt(t(j))))+0.5*exp(x(i)/B^2).*erfc(0.5*(x(i)+t(j))/(B*sqrt(t(j))));

Or if you should insert j+1 in Ur formula do this in line 68

for j=1:M-1 %Time Loop

[Hana Bachi]

I just made the modifications, but the plot I'm getting seems to be totaly wrong as it shows in the attached file

[Abolfazl Chaman Motlagh]

it seems the code is working, but the wrong solution is the result of wrong implementation of Crank-Nicolson method. the U and Ur are not even close to each other in large t (time). in this case i need some time to crack the first numeric part of your code to see what is wrong with it.

[Hana Bachi]

I attached the descritezation I did for crack-Nicolson method. Olease see the attached file

[Hana Bachi]

@Abolfazl Chaman Motlagh

I kept checking the code from the yesterday but I'm still getting the same plots. I think the problem with both the analytical and numerical solution

[Abolfazl Chaman Motlagh]

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]

@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

Answer 2

[Hana Bachi]

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