Fsolve: Issues Solving a system of equations

1 次查看（过去 30 天）
MaxPr
MaxPr 2017-3-6
评论： MaxPr 2017-3-9
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 个评论
Walter Roberson
Walter Roberson 2017-3-8
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
MaxPr 2017-3-9
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!

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

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

请先登录，再进行评论。

Andre
Andre 2017-3-6
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 个评论
MaxPr
MaxPr 2017-3-6
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
Andre 2017-3-6
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by