7 views (last 30 days)

Show older comments

Hey!

So I'm trying to solve this set of non-linear equations. However Fsolve doesn't seem to converge. It goes the right direction however it doesn't seem to quite get there. The solver stops prematurely. I know the math behind the equations is fine, since I already solved it with Mathematica. I just need built it with fsolve.

close all;

%clear all;

eta0=0;

etaInf=9;

deltaEta=0.1; %stepsize eta

N=(etaInf-eta0)/deltaEta+1; %number of nodes +1; later + number of nodes for the belt & the film ;

x0=0;

xInf=1;

deltaX=0.1; %stepsize x

M=(xInf-x0)/deltaX+1; %number of nodes x

anfangswert=horzcat(1.14*ones(M,N),zeros(M,N),0.33*ones(M,N)); %starting values, just zeros for now

%%Aufruf des Solvers

options=optimset('Display','iter','FinDiffType','central','MaxIter',100,'MaxFunEvals',1000000);

Sol = fsolve(@(x)Keller_2D_Diskr_c(x(1:M,1:N),x(1:M,N+1:2*N),...

x(1:M,2*N+1:3*N),N,deltaEta,M,deltaX),anfangswert,options);

%----------------

function fval = Keller_2D_Diskr_c(f,g,H,N,deltaEta,M,deltaX)

% initialising space

fval1=zeros(M,N);fval2=zeros(M,N);fval3=zeros(M,N);

% Define functions of the Form F(X)=0

for i=2:M-1

for j=2:N-1

fval1(i,j)=((f(i,j)-f(i,j-1))/deltaEta)-g(i,j);

fval2(i,j)=((g(i,j)-g(i,j-1))/deltaEta)-H(i,j);

fval3(i,j)=((H(i,j)-H(i,j-1))/deltaEta)+f(i,j)*H(i,j)+...

2*(i/N)*(H(i,j)*((f(i,j)+f(i,j-1)-f(i-1,j)-f(i-1,j-1))...

/2*deltaX)-g(i,j)*((g(i,j)+g(i,j-1)-g(i-1,j)-g(i-1,j-1))...

/2*deltaX));

end

end

%BC's ; favl=[fval1(x,eta),fval2(x,eta),...]

%eta=1.... almost 0

fval(:,1)=[f(:,1);fval2(:,1);1-g(:,1)]; %Boundary values at eta=0, Here f(0)=0, g(0)=1

for j=2:N-1

fval(:,j)=[fval1(:,j-1);-fval2(:,j);fval3(:,j-1)];

end

%eta=N... inf

fval(:,N)=[fval1(:,N-1);g(:,N);fval3(:,N-1)]; %Boundary values at eta=9 g(inf)=0

Any help would be greatly appreciated.

EDIT: So I maybe forgot to mention:

I already tried reducing the stepsizes from deltaEta and deltaX to 0.01. Also I upped the MaxIter and the MaxFunEvals to 1E06. So that shouldn't be it.

My initial guess is actually pretty close, since I'm starting from values from the literature (in this case, the blasius solution of a pulled plate).

Alan Weiss
on 6 Mar 2017

Well, you gave an iteration limit of 100, and that is what stopped the solver. You can pick up the solution process from the final point like this, as recommended by the documentation:

Sol = fsolve(@(x)Keller_2D_Diskr_c(x(1:M,1:N),x(1:M,N+1:2*N),...

x(1:M,2*N+1:3*N),N,deltaEta,M,deltaX),Sol,options);

Or just set a higher value of your iteration limit, maybe 1000 or so.

Alan Weiss

MATLAB mathematical toolbox documentation

Andre
on 6 Mar 2017

Andre
on 6 Mar 2017

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

Start Hunting!