Hi, I am trying to use fsolve for a problem I need to solve.
The code works perfectly for one case (cartesian equations) but when modified to adopt axisymmetric, it does not work.
The error message shows :
"fsolve stopped because the relative size of the current step is less than the
default value of the step size tolerance squared and the vector of function values
is near zero as measured by the default value of the function tolerance."
I have checked the coding several times but I am unable to fix this problem.
I believe it is a problem with optimisation and to do with function solution tolerances. I tried loosening up on MaxFunEvals, MaxIter, which eliminated problem with premature stops, but I still have stalling problem. I also tried TolFun and TolX, unfortunately, no improvement.
If anyone could suggest a solution or point me in the direction, or just anything.....I am struggling big time.
The code of the function is as follows:
% Accepts nx1 vector h, a real x (equidistant mesh size)
% and integer n (# of mesh points)
% nondimensional constant d
% Returns nx1 vector corresponding to A(h)
function [f] = elasticdef(h)
global x X d dt hfirststep ratio_terms gravityterm
n=length(h);
%matrices approximating the derivatives
s0 = ones(n,1); s1 = ones(n1,1);
s2 = ones(n2,1); s3 = ones(n3,1);
D41 = sparse(diag(s3,3) + diag(5*s2,2) + diag(10*s1,1)...
+ diag(10*s0) + diag(5*s1,1) + diag(s2,2));
D42 = sparse(diag(s3,3) + diag(6*s2,2) + diag(15*s1,1)...
+ diag(20*s0) + diag(15*s1,1) + diag(6*s2,2)...
+ diag(s3,3));
D43 = sparse(diag(s2,2) + diag(5*s1,1) + diag(10*s0)...
+ diag(10*s1,1) + diag(5*s2,2) + diag(s3,3));
D21 = sparse(diag(s1,1) + diag(s0));
D22 = sparse(diag(s1,1) + diag(2*s0) + diag(s1,1) );
D23 = sparse(diag(s0) + diag(s1,1));
D1 = sparse(diag(s1,1) + diag(s1,1));
Du = sparse(diag(s1,1)); Dd = sparse(diag(s1,1));
h3 = h.^3;
A4=ratio_terms*((Du*h3).*(D41*h) + (h3).*(D42*h)...
 (Dd*h3).*(D43*h));
A2 = ((Du*h3).*(D21*h) + (h3).*(D22*h)...
 (Dd*h3).*(D23*h));
for i=1:n
A4(i,:)=A4(i,:)./abs(X(i).^3);
A2(i,:)=A2(i,:).*abs(X(i));
end
for i=1:n
A4(i,:) = (A4(i,:)./(2*x^6))./abs(X(i));
A2(i,:) = (A2(i,:)./(2*x^2))./abs(X(i));
end
%%A1 = (D1*h3)/(2*x);
A1 = 0;
%for solving the PDE with 0 inclination
A = A4 + (gravityterm*A2) + A1; %% note the change of sign
A(1) = 0; A(2) = 0;
A(end) = 0; A(end1) = 0;
A(3) = A(3) + ((h3(3)+h3(2) )*h(2))/(2*x^6);
A(end2) = A(end2) + ((h3(end1)+h3(end2))*h(end1))/(2*x^6);
for i=1:n
f(i)=h(i)hfirststep(i)A(i)*dt;
end
================================================
and the code which calls upon the function is as follows:
close all
clear all
N=30;
X=linspace(4,4,N);
h=exp((X.^2));
h=h';
d=1;
x=X(2)X(1);
t = 0; dt = 0.3;
tmax = 2000;
plot(X,h); hold on
global x X d dt hfirststep ratio_terms gravityterm
ratio_terms=1;
hfirststep=h;
gravityterm=1;
for i = 1:4 %two loops are used so the solution can
for j = 1:3 %be plotted only at selected intervals
j , i
t = t+dt;
h=fsolve('elasticdef',hfirststep,optimset('MaxFunEvals',3000000,'MaxIter',400000));
hfirststep=h;
end
plot(X,h); hold on
I=find(h>0.01*max(h));
length_gc(i)=abs(X(max(I))X(min(I)));
end
figure(2)
loglog(length_gc,'k+');
hold on
===============================================
Anyhelp will be greatly appreciated
Thanks a lot for even reading this long post
Thanks a lot
Regards,
Thanks a lot
