what is the fastest appproach to solve a system of non linear equations ?

Hi everyone,

My problem is the following:

i have a rectangular domain discritzed into 1000x100 (MxN) nodal points, where M and N are the total nodal points in the y and x direction respectively.

On each nodal point along the x direction (N=1 to 100) i obtain a non linear system of 1000 algebraic equations which i solve using the function fsolve. so in total for one experimental point i have 100 systems of 1000 equations each. This leads to a simulation time of approximately 15 minutes for just one experimental point.

The issue is that my ultimate goal is to use lsqnonlin to estimate a parameter in my coupled PDE problem by fitting the solution to a total of 75 experimental points. Knowing that just a single simulation of these experimental points would approximately take 19 hours , the fitting procedure will thus take even more time.

So my question is the follwing, is there another function other than fsolve that can be used ? or is there a numerical approach that can minimize my fitting time ?

[Steven Lord]

Can you show how you're setting up your system of equations? Does the function that you're passing into fsolve perform calculations with symbolic variables inside and only substitute values into those equations to generate a numerical result immediately before it returns? If so consider performing the symbolic calculations once and using matlabFunction to generate a function file for fsolve to use that computes the answer solely numerically, no symbolic calculations involved?

Answer 1

[Matt J]

Seems like it would make sense to regard this as a system of 100K simultaneous equations and solve for the PDE parameter with lsqnonlin directly. That way you aren't solving 100 separate times. Also, you seem to imply that you are ultimately pursuing only 1 unknown PDE parameter, so perhaps it would be even better to use fminbnd as opposed to lsqnonlin.

[hamza karim]

Hi both of you,

Thanks for your response and sorry for not elaborating during my problem description.

i attach the system of equations that i am trying to solve. The goal is to estimate the reaction rate constant, kf, and its order, alpha with respect to C1. I was planning to do this by using the lsqnonlin function after i solve my PDEs using an implicit euler scheme.

i cannot consider it as a 100K simultaneous equations because as you will further see in the code, each (n+1)K system depends on the results of the previous one.

An idea that crossed my mind, is to first solve the case of the linear system. This means, that after i obtain my [A , Q] matrices from the Matrix_assembly_BE function, i directly solve for concentration(:,i,j) using the "/" operator. Then once it is computed, after modifyin the fsolve function to take into account only the nonlinear terms and the nodal points that they depend upon, i update the corresponding values in concetration(:,i,j) after using the fsolve function.

Another appraoch is to linearize my source term, and solve the system as a linear system of algebraic equation. However, i am not familiar with the process and any guidance would help me.

The code that i have developped so far is the following.

% data

data=importdata("data.xlsx","data", [2, 75]);

C_initial=rmmissing(data(1:end,1));         % Initial uranium concentration
Keq=rmmissing(data(1:end,2));               % pseudo-equilirbium constant (MODKD correlation)
v=data(1:end,3:4);                          % velocity field 
Cout_exp=data(1:end,5);                     % experimental extraction efficiency (Caq_out/C0)
Daq=rmmissing(data(1:end,6));               % diffusion coefficient in the aqueous phase (m2/s)
NbPoints=rmmissing(data(1:end,7));          % number of experimental points for each concentration
H=205*1E-6;                                 % Channel width
L=1.25*1E-2;                                % Channel Length
Dorg=3*1E-10;                               % Diffusion coefficient in the organic phase
N=100;                                      % Number of nodes in the axial direction
M=1000;                                     % Number of nodes in the y-direction
kf_optimum=1.27E-06;
t=clock;
efficiency=residual_cal_nonlinear(kf_optimum,H,L,N,M,NbPoints,v,Daq,Dorg,Keq,C_initial,Cout_exp);
elapsedtime=etime(clock,t);

This is the main programm, however it is not yet completed where eventually an lsqnonlin function will be used instead. In this example i take the kf value equal to 1.27e-6 and the alpha to be 2. The reason behind this is to just estimate the caclulation time before i do any fitting.

The residual_cal_nolinear will be thus as an input for the lsqnonlin solver. The function calculates the residual between the calculated and experimental points as follows:

function residual=residual_cal_nonlinear(kf,H,L,N,M,NbPoints,v,Daq,Dorg,K,C_initial,Cout_exp)                 

    dy=H/(M-1);
    efficiency=zeros(length(v(1:end,1)),1);
    count=0;
    for n=1:length(C_initial)
    
            
            concentration=zeros(M,N,NbPoints(n));     % Concentration map along the axial direction and y-direction for all experimental points at a specific concentration
                   
            for j=1:NbPoints(n)
                
                C0=zeros(M,1);                              
                C0(1:M/2-1)=C_initial(n);
              
                % solver IMPLICIT EULER
        
                for i=1:N
                    
                    [A,Q]=Matrix_assembly_BE(H,L,N,M,v(count+j,1:end),Daq(n),Dorg,kf,K(n),C0);
                    C_start=ones(M,1);
                    concentration(:,i,j)=fsolve(@nonlinear_system,C_start,[],A,Q,dy,kf,K(n),Daq(n),Dorg);         % solve for the concentration at x+dx
                    C0=concentration(:,i,j);          % re-assigning the initial condition to the last computed solution
                   a=i
                    if(i==N)
                        
                        % Compute extraction efficiency (Caq/Cinitial) at the exit of the channel                    
                        Caq=concentration(1:(M/2-1),N,1);   % Concentration vector in the aqueous phase at the exit of the channel
                        efficiency(count+1)=(2*dy*trapz(Caq))/(H*C_initial(n));
                        
                    end
                end
                                         
            end
            count=count+NbPoints(n);
        
    end
    residual=efficiency-Cout_exp;
end

The main part of the function is the implicit euler solver. Here, the coefficient matrix A and the vector Q are obtained from the Matrix_assembly_BE function. The top (y=0) and bottom (y=H) (see attached file) are obtained using a forward ana backward second order difference schme respectively. Expressions for the internal nodes are obtained by a backward scheme in x and second order central difference scheme in y. The expressions for the interface condition (y=h-) and (y=h+) are obtained using a second order backward and forward scheme respectively.

Note that at the interface i consider the case of alpha =1 (this was the code i orignally developped and i didn't want to modify it).

To take into account the case of the nonlinearity the nonlinear_system function is used where i update the rows that corresponds to the interface condtions by the appropriate expressions.

The functions are as follows:

function [A,Q] = Matrix_assembly_BE(H,L,N,M,v,Daq,Dorg,kf,Keq,C0)

% H: the microchannel width (m)
% L: the microchannel length (m)
% N: Number of nodes in the axial direction
% M: Number of nodes in the y-direction
% v: the velocity vector having a size of 2 v=[vaq vorg]
% Daq,Dorg: Diffusion coefficients
% Keq: The equilibrium constant (Keq=(Keq'x[TBP]^2))

A=zeros(M,M);
Q=zeros(M,1);
dx=L/(N-1);
dy=H/(M-1);      
dy2=dy.^2; 

% Top boundary condition
A(1,1)=3;
A(1,2)=-4;
A(1,3)=1;
Q(1)=0;

% first phase interior nodes
for i=2:M/2-1
    A(i,i)=1/dx+(2*Daq)/(v(1)*dy2);
    A(i,i-1)=-Daq/(v(1)*dy2);
    A(i,i+1)=-Daq/(v(1)*dy2);
    Q(i)=C0(i)/dx;
end

% Uranyl nitrate boundary condition
A(M/2,M/2)=3/(2*dy)+kf/Daq;
A(M/2,M/2-2)=1/(2*dy);
A(M/2,M/2-1)=-2/dy;
A(M/2,M/2+1)=-kf/(Keq*Daq);
Q(M/2)=0;

% Complex boundary condition
A(M/2+1,M/2)=kf/Dorg;
A(M/2+1,M/2+1)=-(3/(2*dy)+kf/(Keq*Dorg));
A(M/2+1,M/2+2)=2/dy;
A(M/2+1,M/2+3)=-1/(2*dy);
Q(M/2+1)=0;

% second phase interior nodes
for i=M/2+2:M-1
    A(i,i)=1/dx+(2*Dorg)/(v(2)*dy2);
    A(i,i-1)=-Dorg/(v(2)*dy2);
    A(i,i+1)=-Dorg/(v(2)*dy2);
    Q(i)=C0(i)/dx;
end

% Bottom boundary condition
A(M,M)=-3;
A(M,M-1)=4;
A(M,M-2)=-1;
Q(M)=0;

end
function F = nonlinear_system(C,A,Q,dy,kf,Keq,Daq,Dorg)
M=length(A(1:end,1));
f=A*C-Q;
f(M/2)=(3/(2*dy))*C(M/2)+(kf/Daq)*C(M/2)^2-(2/dy)*C(M/2-1)+(1/(2*dy))*C(M/2-2)-(kf/(Keq*Daq))*C(M/2+1);
f(M/2+1)=-(3/(2*dy)+kf/(Keq*Dorg))*C(M/2+1)+(kf/Dorg)*C(M/2)^2+(2/dy)*C(M/2+2)-(1/(2*dy))*C(M/2+3);
F=f;
end

[hamza karim]

So an update on the situation, i was able to reduce the computing time from 15 min per experimental point to 37 s by predefining the Jacobian in my objective function. In this case, it would take approximately 5 minutes for a single simulation of all experimental points.

I guess it is reasonable to assume that the fitting procedure would somewhere between 2 and 5 hours. This is reasonable in my case. However if you still have any idea do not hesitate !

[Matt J]

Your follow-up comments haven't really changed my original advice. Since you have a model for the final curve that you're currently fitting, it shouldn't be necessary to solve the ODEs one at a time and post-fit. Instead, it should be possible to develop equations by generating the curve first, based on a hypothetical parameter guess, and then plugging the resulting curve into the ODEs to see if they're satisfied. This would eliminate ODE solving steps altogether.