Crank Nicolson method to solve PDE

Hello, I have the below. When I run it I get the following error message says ''Index in position 1 exceeds array bounds. Index must not exceed 1''

the equation with the boundary and initial conditions are in the attachede file

The code

clear variables
close all
%1. Space steps
xa = 0;
xb = 1;
Lx = 1; %Max value of x
dx =1/40; %space step
N = (xb-xa)/dx; %number of space steps
x=0:dx:Lx;

%2. Time steps
ta = 0;
tb = 0.5;
tf =0.5; %final time
dt =1/3300; %time step
M= (tb-ta)/dt; %number of time of steps
t= 0:dt:tf;
%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


%4.Coefficients of Implicit Method
r=dt/(2*dx^2);
a=1-0.003*r;
b=-0.003*r;
c=2*0.003*r;
e=r*dx/2;
%5. Implicit Method for Au=d Create matriz A and d at level time 1
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); 
A=zeros(N-2,N-2);
d=zeros(N-2,1);
for i=2:N-1 %Space loop
    d(i-1)=(0.003*r)*u(i-1,1)+(0.003*r)*u(i+1,1)+(1-2*0.003*r)*u(i,1)-...
        (r*dx/2)*u(i,1)*u(i+1,1)+(r*dx/2)*u(i,1)*u(i-1,1);
    if i>2 && i<N-1
        A(i-1,i-2:i)=[a b c];
    elseif i==N-1
        A(i-1,end-1:end)=[b c];
        d(i-1)= d(i-1)-a*u(end,2);
    else
        A(i-1,1:2)=[b c];
        d(i-1)= d(i-1)-b*u(1,2);
    end
end
%5. Solving for level time j+1
U=zeros(N,M);
for j=1:M-1 %Time Loop
    u=(A\d)'; %Solving u in the equation Au=d
    U(:,j+1)=[U(1,j+1), u, U(end,j+1)];
    if j~=M-1
        for i=2:N-1 %Space loop
            d(i-1)=(0.003*r)*U(i-1,j+1)+(0.003*r)*U(i+1,j+1)+(1-2*0.003*r)*U(i,j+1)-(r*dx/2)*U(i,j+1)*u(i+1,j+1)+(r*dx/2)*U(i,j+1)*U(i-1,j+1);
            if i==2
                d(i-1) = d(i-1) - b*u(1,j+2);
            elseif i==N-1
                d(i-1) =d(i-1) - a*u(end,j+2);
            end
        end
    end
end
%5. Analytical Solultion and Error
Ur=zeros(N,M);
Error=zeros(N,M);
Exact =@(x,t) phi(x,t);
for j=1:M %Time Loop
    for i=1:N %Space Loop
        
        Error(i,j)=abs(Exact-u(i,j));
    end
end
Exact =@(x,t) phi(x,t);
max_error = max(Error(:,M))
max_error = max(Error(:,round(0.5/dt,0)));
%Plots
nexttile
plot(x,U(:,M));
hold on
plot(x,Ur(:,M));
title('Numerical vs Analytical solution at t=0.85')
xlabel('x-values')
ylabel('U(x,t)')
legend('U @t=0.85 Crank-Nicolson Numerical Solution','U @t=0.85 Analytical Solution')
ylim([0 2.5])
nexttile
plot(x,U(:,round(0.5/dt,0)));
hold on
plot(x,Ur(:,round(0.5/dt,0)));
title('Numerical vs Analytical solution at t=0.5')
xlabel('x-values')
ylabel('U(x,t)')
legend('U @t=0.5 Crank-Nicolson Numerical Solution','U @t=0.5 Analytical Solution')
nexttile
plot(x,Error(:,M));
title('Absolut Error at t=0.85')
xlabel('x')
ylabel('Error')
nexttile
plot(x,Error(:,round(0.5/dt,0)));
title('Absolut error at t=0.5')
xlabel('x')
ylabel('Error')

[KSSV]

Problem is with your variable u. It is of dimension 1x38, i.e. it is an array. But in the loop, where error pops out you are trying to access the element u(3,2), but there is no such element.

[Hana Bachi]

@KSSVHow can I fix that. it says "Error in MathWorkCp4 (line 60)"

Answer 1

[SAI SRUJAN]

Hi Hana,

I understand you are encountering an issue with the Crank-Nicolson method for solving PDEs, specifically an index bounds error within your MATLAB code.

Such errors may arise when an attempt is made to access an array element outside the array's actual bounds for instance, referencing an index that exceeds the array's dimensions.

Upon reviewing the attached code snippet, the error originates from the following line:

d(i-1) = (0.003*r)*U(i-1,j+1) + (0.003*r)*U(i+1,j+1) + (1-2*0.003*r)*U(i,j+1) - (r*dx/2)*U(i,j+1)*u(i+1,j+1) + (r*dx/2)*U(i,j+1)*U(i-1,j+1);

The term `u(i+1,j+1)` seems to be the source of the issue.

Given that `u` is a [1,38] double, attempting to access `u(2,j+1)` when 'i=1' is causing the index bounds error. Please verify the dimensions of the matrix `u` and ensure that any indices referenced do not exceed these dimensions. Adjusting either the indexing or the matrix size should resolve the error.

I hope this helps.