Fsolve: Issues Solving a system of equations

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).

[Walter Roberson]

I did a bit of work on this. You can create symbolic equations from your calculations with a tiny modification to your code:

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

to

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

and then pass in a symbolic array of names to

@(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)

There will be a fair number of 0's in the result, and some duplicates. You can unique the matrix.

After that, you can efficiently solve the first 1636 entries, which are simple linear equations. After that things get messy to solve symbolically, with the computations blowing up quickly, as multiple branches of roots of polynomials need to be included. Possibly constraints would help on that.

This does, though, suggest that you can at least reduce the workload for numeric solving, by solving the easy parts exactly. Reducing from two and a half thousand variables to a thousand variables should help to some degree.

[MaxPr]

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!

Answer 1

[Alan Weiss]

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

[MaxPr]

Well... Thats the last version I added here. I had it run with iteration limit of 1E06 with no luck :/

Answer 2

[Andre]

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.

[MaxPr]

Am I? I was under the impression that the "anfangswert" are my initial values? Because these should be pretty close to the actual solution.

[Andre]

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.