Crank Nicolson method to solve PDE

39 次查看（过去 30 天）
Hana Bachi
Hana Bachi 2022-3-9
回答： SAI SRUJAN 2023-12-26
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')
  2 个评论
KSSV
KSSV 2022-3-9
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
Hana Bachi 2022-3-9
编辑：Hana Bachi 2022-3-9
@KSSVHow can I fix that. it says "Error in MathWorkCp4 (line 60)"

请先登录，再进行评论。

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

回答（1 个）

SAI SRUJAN
SAI SRUJAN 2023-12-26
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.

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by