MATLAB Answers

Fsolve: Issues Solving a system of equations

7 views (last 30 days)
MaxPr
MaxPr on 6 Mar 2017
Commented: MaxPr on 9 Mar 2017
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).
  2 Comments
MaxPr
MaxPr on 9 Mar 2017
Thanks for your suggestion!
I hoped that I missed something in my program, but I do think I have to rethink my way about solving this problem...
Since the part of the program I posted only is a quarter of the whole model.
I'm focusing my efforts now on cycling through a loop with a bvp4c routine. A question might follow ;)
Thanks anyway!

Sign in to comment.

Answers (2)

Alan Weiss
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
  1 Comment
MaxPr
MaxPr on 6 Mar 2017
Well... Thats the last version I added here. I had it run with iteration limit of 1E06 with no luck :/

Sign in to comment.


Andre
Andre on 6 Mar 2017
Can you provide better initial guesses? It seems that you are using zero as initial guess. If you know how the solution looks like, try a guess closer to the solution.
  2 Comments
Andre
Andre on 6 Mar 2017
You have the matrix of zeros, zeros(M,N), inside anfangswert. You can try to initialize with other values. You can also: 1. Try another algorithm for fsolve, using options; 2. Check if you can write the definitions of fval1, 2, 3 in another way; and 3. Use [sol, fval, exitflag] = fsolve to check the message of exitflag and how far you are from the solution with fval.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by