# Error solving Boundary value Problem using bv4pc function

3 views (last 30 days)
Chigga on 12 Jan 2020
Commented: Chigga on 12 Jan 2020
Hi, I am getting the following error while solving this BVP.
Error using bvp4c (line 251)
Unable to solve the collocation equations -- a singular Jacobian encountered.
Error in MaxwellYangHWDiff (line 33)
sol = bvp4c(@(r,y)odefun(r,y,phi,chi),@(ycent,ysurf)bcfun(ycent,ysurf,c1,c2,c3,c4),init);
My code is as follows:
clc
clear all
R = 0.04; % m Radius of the cylinder
f = 2800 * 10^6; % Hertz Frequency of microwave
c = 2.998 * 10^8; % m/s speed of light
mu_o = 1.257 * 10^-6; % Henry/m free space permeability
eps_o = 8.85 * 10^-12; % Frard/m free space permittivity
k1 = 42.6; % dielectric constant
k2 = 13.1; % dielectric loss
P_o = 250; % Watts
a_o = 2*pi*f/c; %free space wavenumber
E_o = sqrt(P_o*pi*a_o*R/(c*eps_o)); % V/m intensity of incident electric field
phi = R^2*(2*pi*f)^2*mu_o*eps_o*k1;
chi = R^2*(2*pi*f)^2*mu_o*eps_o*k2;
nt = 10; %number of nodes
dR = R/(nt-1); %increment
rstar = 0:(dR/R):(R/R);
c1 = R*a_o*((besselj(1,a_o*R)*besselj(0,a_o*R)+bessely(1,a_o*R)*bessely(0,a_o*R))/((besselj(0,a_o*R))^2+(bessely(0,a_o*R))^2));
c2 = 2/(pi*((besselj(0,a_o*R))^2+(bessely(0,a_o*R))^2));
c3 = -(4/pi)*bessely(0,a_o*R)/((besselj(0,a_o*R))^2+(bessely(0,a_o*R))^2);
c4 = -(4/pi)*besselj(0,a_o*R)/((besselj(0,a_o*R))^2+(bessely(0,a_o*R))^2);
init = bvpinit(linspace(0,0.04,10),[0 0 0 0]);
sol = bvp4c(@(r,y)odefun(r,y,phi,chi),@(ycent,ysurf)bcfun(ycent,ysurf,c1,c2,c3,c4),init);
x = linspace(0,0.04,100); BS=deval(sol,x);
plot(x,BS(1,:))
function dydr = odefun(r,y,phi,chi)
dydr=zeros(4,1);
dydr(2)= -((1/r)*y(2)+phi*y(1)-chi*y(3));
dydr(4)= -((1/r)*y(4)+phi*y(3)-chi*y(1));
end
function res = bcfun(ycent,ysurf,c1,c2,c3,c4)
res =[ycent(2); ycent(4); ysurf(2)+c1*ysurf(1)+c2*ysurf(3)-c3; ysurf(4)+c1*ysurf(3)-c2*ysurf(1)-c4];
end

darova on 12 Jan 2020
Try r = 1e-3
I think you forgot about dydr(1) and dydr(3)
function dydr = odefun(r,y,phi,chi)
dydr=zeros(4,1);
% dydr(1) = y(2);
dydr(2)= -((1/r)*y(2)+phi*y(1)-chi*y(3));
% dydr(3) = y(4);
dydr(4)= -((1/r)*y(4)+phi*y(3)-chi*y(1));
end

Chigga on 12 Jan 2020
hi, thank you for your answer. I made the changes you suggested. The error I mentioned earlier no longer exists. However, I am getting the following error:
Attempted to access x(1); index out of bounds because numel(x)=0.
Error in bvpinit (line 71)
if x(1) == x(N)
Error in MaxwellYangHWDiff (line 28)
init = bvpinit(linspace(0,0.04,1e-3),[0 0 0 0]);
I tried using different intial points, but it still gives me this error.
Also, there is a correction.
function dydr = odefun(r,y,phi,chi)
dydr=zeros(4,1);
dydr(1) = y(2);
dydr(2)= -((1/r)*y(2)+phi*y(1)-chi*y(3));
dydr(3) = y(4);
dydr(4)= -((1/r)*y(4)+phi*y(3)+chi*y(1));
end
darova on 12 Jan 2020
I tried this
init = bvpinit(linspace(eps,0.4,100),[0 0 0 0]);
sol = bvp4c(@odefun, @bcfun, init);
plot(sol.x,sol.y)
And get this
See attached script
Chigga on 12 Jan 2020
Thank you! I really appreciate your help!

Translated by