201 views (last 30 days)

Show older comments

I am solving the laplace equation in spherical coordinates with a cell (spherical capacitor) at the origin. The equation becomes,

After discretizing using the finite difference method and applying the appropriate boundary conditions, I am getting the correct answer off by ~ 0.3% , I am also recieving the warning message that the Matrix is close to singular. Is RCOND = 9.735910e-18 small enough that I can ignore this warning? How does this warning effect my solution?

close all

clear all

clc

%% setting up A matrix

a = 50e-6; % cell radius

dth = pi/128; % angle step size

dp = 3*a/60; % radius step size

dt = 1e-8; % time step

angles = dth/2:dth:2*pi-dth/2; % angle values

radii = 0:dp:3*a; % radii values

E = 40e3 ; % applid field strength

C = 1e-2; % Capacitance

g = 2; % Conductance

si = 0.455; % intracellular conductivity

se = 5; %extracellular conductivity

Vrest = -0.08;

% constructing 'A' Matrix

A = zeros(15616,15872);

k = 256;

for i = 1:60

for j = 1:256

if j == 1 && i ~= 21 && i~=20 && i~=60

A(k*(i-1)+j,k*i+j) = -2/dp^2 - 2/(((i*dp)^2)*(dth^2)); %U(i,j)

A(k*(i-1)+j,i*k+1+j) = 1/(((i*dp)^2)*(dth^2)) + cot((2*pi/256)*(j-1+0.5))/(((i*dp)^2)*2*dth) ; %U(i,j+1)

A(k*(i-1)+j,(i+1)*k) = 1/(((i*dp)^2)*(dth^2)) - cot((2*pi/256)*(j-1+0.5))/(((i*dp)^2)*2*dth); %U(i,j-1)

A(k*(i-1)+j,(i-1)*k + j) = 1/dp^2 - 2/(2*i*dp*dp); %U(i-1,j)

A(k*(i-1)+j,k*(i+1)+j) = 1/dp^2 + 2/(2*i*dp*dp); %U(i+1,j)

elseif j == k && i ~= 21 && i~=20 && i~=60

A(k*(i-1)+j,k*i+j) = -2/dp^2 - 2/(((i*dp)^2)*(dth^2));

A(k*(i-1)+j,i*k+1) = 1/(((i*dp)^2)*(dth^2)) + cot((2*pi/256)*(j-1+0.5))/(((i*dp)^2)*2*dth) ;

A(k*(i-1)+j,i*k-1+j) = 1/(((i*dp)^2)*(dth^2)) - cot((2*pi/256)*(j-1+0.5))/(((i*dp)^2)*2*dth);

A(k*(i-1)+j,(i-1)*k+j) = 1/dp^2 - 2/(2*i*dp*dp);

A(k*(i-1)+j,k*(i+1)+j) = 1/dp^2 + 2/(2*i*dp*dp);

elseif i ~= 21 && i~=20 && i~=60

A(k*(i-1)+j,k*i+j) = -2/dp^2 - 2/(((i*dp)^2)*(dth^2));

A(k*(i-1)+j,i*k+1+j) = 1/(((i*dp)^2)*(dth^2)) + cot((2*pi/256)*(j-1+0.5))/(((i*dp)^2)*2*dth) ;

A(k*(i-1)+j,i*k-1+j) = 1/(((i*dp)^2)*(dth^2)) - cot((2*pi/256)*(j-1+0.5))/(((i*dp)^2)*2*dth);

A(k*(i-1)+j,(i-1)*k+j) = 1/dp^2 - 2/(2*i*dp*dp);

A(k*(i-1)+j,k*(i+1)+j) = 1/dp^2 + 2/(2*i*dp*dp);

elseif i == 21

A(k*(i-1)+j,k*i+j) = 3*se/(2*dp); %U(i,j)

A(k*(i-1)+j,k*(i+1)+j) = -4*se/(2*dp); %U(i+1,j)

A(k*(i-1)+j,k*(i+2)+j) = se/(2*dp) ; %U(i+2,j)

A(k*(i-1)+j,15616+j) = -C/dt -g; %U(Vm,j)

elseif i == 20

A(k*(i-1)+j,k*i+j) = -3*si/(2*dp); %U(i,j)

A(k*(i-1)+j,k*(i-1)+j) = 4*si/(2*dp); %U(i-1,j)

A(k*(i-1)+j,k*(i-2)+j) = -si/(2*dp) ; %U(i-2,j)

A(k*(i-1)+j,15616+j) = -C/dt -g; %U(Vm,j)

elseif i == 60

A(k*(i)+j,k*(i+1)+j) = 1; %U(i,j)

A(k*(i)+j,k*20+j) = -1; %U(20,j)

A(k*(i)+j,k*21+j) = 1; %U(21,j)

end

end

end

% origin is average of surounding disk

A = [zeros(256,15872) ; A ] ;

A(1:256,1:256) = diag(-256*ones(1,256));

A(1:256,257:512) = ones(256,256);

A(15361:end-256,15361:end-256) = diag(ones(1,256)); % Boundary Condition E_app

b = zeros(15872,1);

b(15361:end-256) = -E*3*a*cos(angles); % Boundary Condition E_app

b(5121:5376) = -1*(Vrest*(C/dt) + g*Vrest);

b(5377:5632) = -1*(Vrest*(C/dt) + g*Vrest);

dA = decomposition(A);

x = dA\b;

X = transpose(reshape(x,[256,62]));

Vm = zeros(1000,256);

Vm(1,:) = -0.08;

for t = 1:1000

Vm(t+1,:) = X(62,:);

b(5121:5376) = -1*((transpose(X(62,:)))*(C/dt) + g*Vrest);

b(5377:5632) = -1*((transpose(X(62,:)))*(C/dt) + g*Vrest);

x = dA\b;

X = transpose(reshape(x,[256,62]));

end

imagesc(X(1:end-1,:))

Christine Tobler
on 8 Jan 2021

Edited: Christine Tobler
on 8 Jan 2021

I tried plotting the following two things:

>> figure; semilogy(max(abs(A), [], 1))

>> figure; semilogy(max(abs(A), [], 2))

From these, it looks like the first 200 columns and the last few hundreds rows of A have much smaller scale than the others. This will usually lead to a bad condition number - if a row or column of A was exactly zero, that matrix A would be singular.

If it makes sense for your problem to change the scaling on those rows and/or columns, that is likely to improve the conditioning of A. Possibly changing the scaling of the boundary conditions would solve the problem.

Keep in mind this isn't guaranteed: If a matrix has some rows or columns with much smaller scale (about 1e-15) than others, this means it has a bad condition number - but if all rows and columns have similar scaling, that's not enough to make it well conditioned in general (think of the matrix of all ones for example, which is singular).

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

Start Hunting!