Singular Jacobian Encountered Error- bvp4c
1 view (last 30 days)
Show older comments
function sol = singular_flameballT1
format long e clc clear all global beta global gamma global Le global sl global S len1 = 0.0; % solving as a singular problem based on bvpset reference from matlab site len2 = 120;
beta = 10; gamma = 5; % gamma = q Le = 1.0; %DO NOT MODIFY THIS VALUE,CHANGE BELOW 4 lines sl = 0.977865;
x = len1:0.01:len2;
% generating the constant matrix for solving as a singular problem S = [0 0 0 -2.0];
options = bvpset('NMax',1200001,'SingularTerm',S,'FJacobian',@fjacob,'BCJacobian',@bcjacob,'Vectorized','on');
%load radLe1beta10len120step.mat;
solinit=bvpinit(x,@mat4init); sol=bvp4c(@eigen_value_sol,@bc,solinit,options); fprintf('First sol done!\n');
output(sol,len1,len2);
end
function [] = output(sol,len1,len2) global gamma global beta global Le global sl
figure(1) plot(sol.x,sol.y(1,:),'.');
hold on
Temp(1,:) = sol.y(1,:); Y(1,:) = 1-sol.y(1,:); plot(sol.x,Y(1,:),'r'); % calculating omega [~,b]=size(Temp); hold on m=zeros(1,b); for i=1:b m(1,i)=((1+gamma)/(1+(gamma*Temp(1,i)/Le)))*Y(1,i)*((beta^2)/(2*(sl^2)*Le))*exp((beta*(Temp(1,i)/Le-1))/((1+gamma*Temp(1,i)/Le)/(gamma+1))); end
fprintf('The omega at 0 is %13.9f.\n',m(1,1));
hold on %figure(3) plot(sol.x,m,'g'); % axis([len1 len2 -0.2 2]) legend('T','Y','w','location','NorthEastOutside'); xlabel('x') ylabel('solution') hold on
end
function f = eigen_value_sol(x,y)
global beta global gamma global Le global sl f=[ y(2,:) -((1+gamma)./(1+((gamma/Le)*y(1,:)))).*(1-y(1,:)).*(((beta^2)/(2*(sl^2)))*exp((beta*(y(1,:)/Le-1))./((1+(gamma/Le)*y(1,:))/(gamma+1)))) ]; end
function fjac = fjacob(x,y) global gamma global beta global sl global Le t1 = -((beta^2)/(2*(sl^2)))*(1+gamma); t3 = exp((beta*(y(1)/Le-1))/((1+gamma*y(1)/Le)/(gamma+1))); t2 = (1-y(1))/(1+(gamma*y(1)/Le)); dev1 = t1*t3*(1/(1+(gamma*y(1)/Le))^2)*(-1-gamma/Le) + t1*t2*t3*beta*(1+gamma)*(1/(1+(gamma*y(1)/Le))^2)*(1/Le)*(1+gamma);
fjac = [0 1 dev1 0];
end
function [bcajac,bcbjac] = bcjacob(x,y) % global gamma % global beta % global sl % global Le
bcajac = [0 1 0 0];
bcbjac = [0 0 1 0];
end
%boundary conditions function res=bc(ya,yb) %global Le %Dont change res=[ya(2) yb(1)];
end
%starting solution function yinit = mat4init(x) global Le c=4.0; a=2.0; % this is where the intersection happpens
% Le scaling and without derivatives yinit = [((1+erf(c*(a-x)))/2) 0];
% yinit = [((1+erf(c*(a-x)))/2) % ((-c/(pi^(1/2))*exp(-(c*(a-x))^2))) % ]; end
Any suggestions on why the singular jacobian error might occur in above code? Thanks in advance!
0 Comments
Answers (0)
See Also
Categories
Find more on Create Model Web Views in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!