MATLAB Central Newsreader - Solving implicit solution

Solving implicit solution

Eng M — Sun, 05 Feb 2012 15:07:10 +0000

I need to find d for the following code and i am getting an error, what is the problem?

fy=40;
cover=1.25;
ql=1;
for k=1:11
    Lx(k)=120+(k-1)*24;
for l=1:11
    Ly(l)=120+(l-1)*24;
for i=1:11
    Cx(i)=8+(i-1)*4;
for j=1:11
    Cy(j)=8+(j-1)*4;

if (Cx(i)>=Cy(j))
    B(i,j)=Cx(i)/Cy(j);
else B(i,j)=Cy(j)/Cx(i);
end
d=zeros([11 11 11 11]);
A(i,j,k,l)=(Lx(k)+Ly(l))/(Cx(i)+Cy(j));
G1(i,j,k,l)=d(i,j,k,l)*(Cx(i)+Cy(j)+2*d(i,j,k,l))/((Lx(k)*Ly(l)-(Cx(i)+d(i,j,k,l))*(Cy(j)+d(i,j,k,l)))*(1+(3/2)*sqrt((Cy(j)+d(i,j,k,l))/(Cx(i)+d(i,j,k,l)))))*0.336*Ly(l)*(Lx(k)-Cx(i))^2/(d(i,j,k,l)^3+(Cx(i)+d(i,j,k,l))^2*d(i,j,k,l)+3*d(i,j,k,l)*(Cy(j)+d(i,j,k,l))*(Cx(i)+d(i,j,k,l)));
G2(i,j,k,l)=d(i,j,k,l)*(Cx(i)+Cy(j)+2*d(i,j,k,l))/((Lx(k)*Ly(l)-(Cx(i)+d(i,j,k,l))*(Cy(j)+d(i,j,k,l)))*(1+(3/2)*sqrt((Cx(i)+d(i,j,k,l))/(Cy(j)+d(i,j,k,l)))))*0.336*Lx(k)*(Ly(l)-Cy(j))^2/(d(i,j,k,l)^3+(Cy(j)+d(i,j,k,l))^2*d(i,j,k,l)+3*d(i,j,k,l)*(Cx(i)+d(i,j,k,l))*(Cy(j)+d(i,j,k,l)));
delta1(i,j,k,l)= 1+((5/32)*(G1(i,j,k,l)+G2(i,j,k,l)));

if ((2+4/B(i,j))<= ((40*d(i,j,k,l))/((2*(Cx(i)+d(i,j,k,l))+2*(Cy(j)+d(i,j,k,l))))+2)) && ((2+4/B(i,j))<=4)
    Vc(i,j,k,l)=(2+4/B(i,j));
elseif (((40*d(i,j,k,l))/((2*(Cx(i)+d(i,j,k,l))+2*(Cy(j)+d(i,j,k,l))))+2)<= (2+4/B(i,j))) && (((40*d(i,j,k,l))/((2*(Cx(i)+d(i,j,k,l))+2*(Cy(j)+d(i,j,k,l))))+2)<=4)
    Vc(i,j,k,l)= ((40*d(i,j,k,l))/((2*(Cx(i)+d(i,j,k,l))+2*(Cy(j)+d(i,j,k,l))))+2);
elseif (4<= ((40*d(i,j,k,l))/((2*(Cx(i)+d(i,j,k,l))+2*(Cy(j)+d(i,j,k,l))))+2)) && (4<=(2+4/B(i,j)))
    Vc(i,j,k,l)= 4;

end
solve('delta1*6.4*ql*(Lx(k)*Ly(l)-(Cx(i)+d(i,j,k,l))*(Cy(j)+d(i,j,k,l)))/(d(i,j,k,l)*(2*Cx(i)+2*Cy(j)+4*d(i,j,k,l)))=Vc','d(i,j,k,l)');

end
end
end
end
end
a=reshape (A,14641,1);
D=reshape (d,14641,1);

Thank you in advance

Re: Solving implicit solution

Roger Stafford — Sun, 05 Feb 2012 19:34:10 +0000

"Eng M" wrote in message <jgm5ut$o0q$1@newscl01ah.mathworks.com>...
> I need to find d for the following code and i am getting an error, what is the problem?
> solve('delta1*6.4*ql*(Lx(k)*Ly(l)-(Cx(i)+d(i,j,k,l))*(Cy(j)+d(i,j,k,l)))/(d(i,j,k,l)*(2*Cx(i)+2*Cy(j)+4*d(i,j,k,l)))=Vc','d(i,j,k,l)');
- - - - - - - - - - - - -
That's just a quadratic equation in d. Its coefficients are complicated but there is a simple formula for solving it learned by high school students. Why are you using the sophisticated 'solve' for such an elementary purpose?

Roger Stafford

Re: Solving implicit solution

Eng M — Mon, 06 Feb 2012 14:46:10 +0000

"Roger Stafford" wrote in message <jgmlji$98v$1@newscl01ah.mathworks.com>...
> "Eng M" wrote in message <jgm5ut$o0q$1@newscl01ah.mathworks.com>...
> > I need to find d for the following code and i am getting an error, what is the problem?
> > solve('delta1*6.4*ql*(Lx(k)*Ly(l)-(Cx(i)+d(i,j,k,l))*(Cy(j)+d(i,j,k,l)))/(d(i,j,k,l)*(2*Cx(i)+2*Cy(j)+4*d(i,j,k,l)))=Vc','d(i,j,k,l)');
> - - - - - - - - - - - - -
> That's just a quadratic equation in d. Its coefficients are complicated but there is a simple formula for solving it learned by high school students. Why are you using the sophisticated 'solve' for such an elementary purpose?
>
> Roger Stafford
--------------------------------------
It is not quadratic equation where delta1 is also a function of d!!! and as you can see i have a function of d=minimum between 3 functions one of them also is a function of d !!!

Re: Solving implicit solution

Roger Stafford — Mon, 06 Feb 2012 17:29:09 +0000

"Eng M" wrote in message <jgop3i$nkq$1@newscl01ah.mathworks.com>...
> It is not quadratic equation where delta1 is also a function of d!!! and as you can see i have a function of d=minimum between 3 functions one of them also is a function of d !!!
- - - - - - - - -
  As you have things set up, 'solve' thinks it is dealing with the equivalent of a quadratic equation. It does not know about your apparent intent that G1, G2, delta1, and Vc should all vary as d varies. In your code you have set d to zero in computing them. These are the values 'solve' uses in attempting to derive a solution, not the variable values.

  It would be extremely difficult to express your problem in a form that would be suitable for 'solve', given the complexity I see here. You should be using 'fzero' where provision is made for the functionality you have in mind.

  My apologies for my earlier, overly hasty analysis of your problem.

Roger Stafford

Re: Solving implicit solution

Eng M — Mon, 06 Feb 2012 18:22:09 +0000

"Roger Stafford" wrote in message <jgp2l5$spg$1@newscl01ah.mathworks.com>...
> "Eng M" wrote in message <jgop3i$nkq$1@newscl01ah.mathworks.com>...
> > It is not quadratic equation where delta1 is also a function of d!!! and as you can see i have a function of d=minimum between 3 functions one of them also is a function of d !!!
> - - - - - - - - -
> As you have things set up, 'solve' thinks it is dealing with the equivalent of a quadratic equation. It does not know about your apparent intent that G1, G2, delta1, and Vc should all vary as d varies. In your code you have set d to zero in computing them. These are the values 'solve' uses in attempting to derive a solution, not the variable values.
>
> It would be extremely difficult to express your problem in a form that would be suitable for 'solve', given the complexity I see here. You should be using 'fzero' where provision is made for the functionality you have in mind.
>
> My apologies for my earlier, overly hasty analysis of your problem.
>
> Roger Stafford
----------------------
Thank you very much I tried to use fsolve but i am still getting errors:
This is my new code i tried to combine all G1,G2, and delta 1 in one equation:

fy=40;
cover=1.25;
ql=1;
for k=1:11
    Lx(k)=120+(k-1)*24;
for l=1:11
    Ly(l)=120+(l-1)*24;
for i=1:11
    Cx(i)=8+(i-1)*4;
for j=1:11
    Cy(j)=8+(j-1)*4;

if (Cx(i)>=Cy(j))
    B(i,j)=Cx(i)/Cy(j);
else B(i,j)=Cy(j)/Cx(i);
end
A(i,j,k,l)=(Lx(k)+Ly(l))/(Cx(i)+Cy(j));

fsolve(@(d)((1+((5/32)*((d*(Cx(i)+Cy(j)+2*d)/((Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))*(1+(3/2)*sqrt((Cy(j)+d)/(Cx(i)+d))))*0.336*Ly(l)*(Lx(k)-Cx(i))^2/(d^3+(Cx(i)+d)^2*d+3*d*(Cy(j)+d)*(Cx(i)+d)))+(d*(Cx(i)+Cy(j)+2*d)/((Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))*(1+(3/2)*sqrt((Cx(i)+d)/(Cy(j)+d))))*0.336*Lx(k)*(Ly(l)-Cy(j))^2/(d^3+(Cy(j)+d)^2*d+3*d*(Cx(i)+d)*(Cy(j)+d))))))*6.4*ql*(Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))/(d*(2*Cx(i)+2*Cy(j)+4*d))-min([2+4/B(i,j);(40*d)/((2*(Cx(i)+d)+2*(Cy(j)+d)))+2;4])),3);
D=d(i,j,k,l);
end
end
end
end

a=reshape (A,14641,1);
depth=reshape(D,14641,1);

Could you please help me in this?
Thank you in advance

Re: Solving implicit solution

Roger Stafford — Mon, 06 Feb 2012 19:43:11 +0000

"Eng M" wrote in message <jgp5oh$am9$1@newscl01ah.mathworks.com>...
> Thank you very much I tried to use fsolve but i am still getting errors:
> This is my new code i tried to combine all G1,G2, and delta 1 in one equation:
>
> fy=40;
> cover=1.25;
> ql=1;
> for k=1:11
> Lx(k)=120+(k-1)*24;
> for l=1:11
> Ly(l)=120+(l-1)*24;
> for i=1:11
> Cx(i)=8+(i-1)*4;
> for j=1:11
> Cy(j)=8+(j-1)*4;
>
> if (Cx(i)>=Cy(j))
> B(i,j)=Cx(i)/Cy(j);
> else B(i,j)=Cy(j)/Cx(i);
> end
> A(i,j,k,l)=(Lx(k)+Ly(l))/(Cx(i)+Cy(j));
>
> fsolve(@(d)((1+((5/32)*((d*(Cx(i)+Cy(j)+2*d)/((Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))*(1+(3/2)*sqrt((Cy(j)+d)/(Cx(i)+d))))*0.336*Ly(l)*(Lx(k)-Cx(i))^2/(d^3+(Cx(i)+d)^2*d+3*d*(Cy(j)+d)*(Cx(i)+d)))+(d*(Cx(i)+Cy(j)+2*d)/((Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))*(1+(3/2)*sqrt((Cx(i)+d)/(Cy(j)+d))))*0.336*Lx(k)*(Ly(l)-Cy(j))^2/(d^3+(Cy(j)+d)^2*d+3*d*(Cx(i)+d)*(Cy(j)+d))))))*6.4*ql*(Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))/(d*(2*Cx(i)+2*Cy(j)+4*d))-min([2+4/B(i,j);(40*d)/((2*(Cx(i)+d)+2*(Cy(j)+d)))+2;4])),3);
> D=d(i,j,k,l);
> end
> end
> end
> end
>
> a=reshape (A,14641,1);
> depth=reshape(D,14641,1);
>
> Could you please help me in this?
> Thank you in advance
- - - - - - - - -
In the assignment "D=d(i,j,k,l)" it looks as if you have things backwards. The quantity 'fsolve' (hopefully) solves will be a simple scalar quantity and matlab would not know what to do with something like d(i,j,k,l). I think you should be writing "d(i,j,k,l)=D" with D being the variable solved for by 'solve'.

A piece of advice: When presenting a problem of this complex nature to this newsreader audience, you should reduce it to its simplest form. Forget the fact that this is occurring within a set of nested for-loops, which has nothing to do with your difficulties, enabling you to leave out all the indices such as i in Cx(i) and just concentrate on the problem of solving a single implicit equation. Try lumping the various known values together into single quantities to get as simple an expression as possible. You could have left out altogether the computation of A which has nothing to do with your problem. Also it is not necessary that we know the details of how Cx, Cy, and B are computed, but simply the fact that they are known. Make it clear to us (and to 'solve') what your knowns and unknowns are to be. Your original formulation of the input to 'solve' contained over a thousand
characters and the more recent one is not much simpler, which is certainly discouraging to volunteers like myself with limited time and energy.

Roger Stafford

Re: Solving implicit solution

Eng M — Mon, 06 Feb 2012 20:40:10 +0000

"Roger Stafford" wrote in message <jgpagf$r51$1@newscl01ah.mathworks.com>...
> "Eng M" wrote in message <jgp5oh$am9$1@newscl01ah.mathworks.com>...
> > Thank you very much I tried to use fsolve but i am still getting errors:
> > This is my new code i tried to combine all G1,G2, and delta 1 in one equation:
> >
> > fy=40;
> > cover=1.25;
> > ql=1;
> > for k=1:11
> > Lx(k)=120+(k-1)*24;
> > for l=1:11
> > Ly(l)=120+(l-1)*24;
> > for i=1:11
> > Cx(i)=8+(i-1)*4;
> > for j=1:11
> > Cy(j)=8+(j-1)*4;
> >
> > if (Cx(i)>=Cy(j))
> > B(i,j)=Cx(i)/Cy(j);
> > else B(i,j)=Cy(j)/Cx(i);
> > end
> > A(i,j,k,l)=(Lx(k)+Ly(l))/(Cx(i)+Cy(j));
> >
> > fsolve(@(d)((1+((5/32)*((d*(Cx(i)+Cy(j)+2*d)/((Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))*(1+(3/2)*sqrt((Cy(j)+d)/(Cx(i)+d))))*0.336*Ly(l)*(Lx(k)-Cx(i))^2/(d^3+(Cx(i)+d)^2*d+3*d*(Cy(j)+d)*(Cx(i)+d)))+(d*(Cx(i)+Cy(j)+2*d)/((Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))*(1+(3/2)*sqrt((Cx(i)+d)/(Cy(j)+d))))*0.336*Lx(k)*(Ly(l)-Cy(j))^2/(d^3+(Cy(j)+d)^2*d+3*d*(Cx(i)+d)*(Cy(j)+d))))))*6.4*ql*(Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))/(d*(2*Cx(i)+2*Cy(j)+4*d))-min([2+4/B(i,j);(40*d)/((2*(Cx(i)+d)+2*(Cy(j)+d)))+2;4])),3);
> > D=d(i,j,k,l);
> > end
> > end
> > end
> > end
> >
> > a=reshape (A,14641,1);
> > depth=reshape(D,14641,1);
> >
> > Could you please help me in this?
> > Thank you in advance
> - - - - - - - - -
> In the assignment "D=d(i,j,k,l)" it looks as if you have things backwards. The quantity 'fsolve' (hopefully) solves will be a simple scalar quantity and matlab would not know what to do with something like d(i,j,k,l). I think you should be writing "d(i,j,k,l)=D" with D being the variable solved for by 'solve'.
>
> A piece of advice: When presenting a problem of this complex nature to this newsreader audience, you should reduce it to its simplest form. Forget the fact that this is occurring within a set of nested for-loops, which has nothing to do with your difficulties, enabling you to leave out all the indices such as i in Cx(i) and just concentrate on the problem of solving a single implicit equation. Try lumping the various known values together into single quantities to get as simple an expression as possible. You could have left out altogether the computation of A which has nothing to do with your problem. Also it is not necessary that we know the details of how Cx, Cy, and B are computed, but simply the fact that they are known. Make it clear to us (and to 'solve') what your knowns and unknowns are to be. Your original formulation of the input to 'solve' contained over a thousand
> characters and the more recent one is not much simpler, which is certainly discouraging to volunteers like myself with limited time and energy.
>
> Roger Stafford
--------
Thank you very much sir, and sorry for disturbance. This is the first time that i use this newsreader, and i just thought writing the whole code will make it easier to understand what i wanted excatly (to show that i want to find solution for function inside 4 loops). Thank you for your time and efford.

Re: Solving implicit solution

Roger Stafford — Tue, 07 Feb 2012 01:42:10 +0000

"Eng M" wrote in message <jgpdra$9eg$1@newscl01ah.mathworks.com>...
> Thank you very much sir, and sorry for disturbance. This is the first time that i use this newsreader, and i just thought writing the whole code will make it easier to understand what i wanted excatly (to show that i want to find solution for function inside 4 loops). Thank you for your time and efford.
- - - - - - - - -
I should point out that since both 'fsolve' and 'fzero' require an initial estimate be entered as an input, perhaps you can take advantage of one result to serve as a smart estimate for the next call. This assumes that your Lx, Ly, Cx, and Cy values will not usually change very much from one time to another and therefore hopefully the solution would not change greatly. This might save on execution time by avoiding wide swings in the solution search which these routines must perform.

But of course before such a refinement you must first manage to get some one single solution.

Roger Stafford

Re: Solving implicit solution

Eng M — Tue, 07 Feb 2012 14:43:09 +0000

"Roger Stafford" wrote in message <jgpvhi$4j9$1@newscl01ah.mathworks.com>...
> "Eng M" wrote in message <jgpdra$9eg$1@newscl01ah.mathworks.com>...
> > Thank you very much sir, and sorry for disturbance. This is the first time that i use this newsreader, and i just thought writing the whole code will make it easier to understand what i wanted excatly (to show that i want to find solution for function inside 4 loops). Thank you for your time and efford.
> - - - - - - - - -
> I should point out that since both 'fsolve' and 'fzero' require an initial estimate be entered as an input, perhaps you can take advantage of one result to serve as a smart estimate for the next call. This assumes that your Lx, Ly, Cx, and Cy values will not usually change very much from one time to another and therefore hopefully the solution would not change greatly. This might save on execution time by avoiding wide swings in the solution search which these routines must perform.
>
> But of course before such a refinement you must first manage to get some one single solution.
>
> Roger Stafford
---------------------
I am trying to define d as an array by many methods but still i am getting an error, it seems that to use 'fsolve' d should not be an array but a single value. So, what should I do to let it work in a loop?

Re: Solving implicit solution

Roger Stafford — Tue, 07 Feb 2012 17:35:10 +0000

"Eng M" wrote in message <jgrd9t$q68$1@newscl01ah.mathworks.com>...
> I am trying to define d as an array by many methods but still i am getting an error, it seems that to use 'fsolve' d should not be an array but a single value. So, what should I do to let it work in a loop?
- - - - - - - - - - -
Yes you will need a loop eventually. However, before bothering to write such a loop you need to concentrate on getting your code working for a single set of Cx, etc. values. From what I can see you have not succeeded in that yet. Using one set of Lx, Ly, Cx, and Cy values you need to coerce your code into producing a valid solution. After that, putting it all in a loop is easy.

Roger Stafford

Re: Solving implicit solution

Eng M — Tue, 07 Feb 2012 17:54:10 +0000

"Roger Stafford" wrote in message <jgrnce$2s6$1@newscl01ah.mathworks.com>...
> "Eng M" wrote in message <jgrd9t$q68$1@newscl01ah.mathworks.com>...
> > I am trying to define d as an array by many methods but still i am getting an error, it seems that to use 'fsolve' d should not be an array but a single value. So, what should I do to let it work in a loop?
> - - - - - - - - - - -
> Yes you will need a loop eventually. However, before bothering to write such a loop you need to concentrate on getting your code working for a single set of Cx, etc. values. From what I can see you have not succeeded in that yet. Using one set of Lx, Ly, Cx, and Cy values you need to coerce your code into producing a valid solution. After that, putting it all in a loop is easy.
>
> Roger Stafford
-------------
It is working for one set of Cx,Cy,Lx, and Ly, but it is not working when it comes to full set because fsolve is not working for multi-variables but only one set of Cx,Cy,Lx and Ly as you said!

Re: Solving implicit solution

Eng M — Tue, 07 Feb 2012 18:14:09 +0000

"Eng M" wrote in message <jgrog2$6jt$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <jgrnce$2s6$1@newscl01ah.mathworks.com>...
> > "Eng M" wrote in message <jgrd9t$q68$1@newscl01ah.mathworks.com>...
> > > I am trying to define d as an array by many methods but still i am getting an error, it seems that to use 'fsolve' d should not be an array but a single value. So, what should I do to let it work in a loop?
> > - - - - - - - - - - -
> > Yes you will need a loop eventually. However, before bothering to write such a loop you need to concentrate on getting your code working for a single set of Cx, etc. values. From what I can see you have not succeeded in that yet. Using one set of Lx, Ly, Cx, and Cy values you need to coerce your code into producing a valid solution. After that, putting it all in a loop is easy.
> >
> > Roger Stafford
> -------------
> It is working for one set of Cx,Cy,Lx, and Ly, but it is not working when it comes to full set because fsolve is not working for multi-variables but only one set of Cx,Cy,Lx and Ly as you said!
----------------------
For single set i am getting the following:

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

Re: Solving implicit solution

Christopher Creutzig — Mon, 20 Feb 2012 12:28:40 +0000

On 05.02.12 16:07, Eng M wrote:
> I need to find d for the following code and i am getting an error, what is the problem?

> solve('delta1*6.4*ql*(Lx(k)*Ly(l)-(Cx(i)+d(i,j,k,l))*(Cy(j)+d(i,j,k,l)))/(d(i,j,k,l)*(2*Cx(i)+2*Cy(j)+4*d(i,j,k,l)))=Vc','d(i,j,k,l)');

In addition to a number of things Roger already said: You are using
strings in the input to solve. That means your values for Lx, Ly, etc.
are not used.

Oh, and you didn't tell us which error you got. Typically, the error
message contains at least an attempt at something helpful and they are
often worth reading carefully.

Christopher

Re: Solving implicit solution

Marc — Tue, 21 Feb 2012 02:16:10 +0000

You were close, just a couple of formatting changes seem to get me an answer.

I made the following changes to your previous code....

In the way beginning, before ANY looping, i would preallocate all of the vectors and arrays. Also, set oldD = 3.;

For fsolve, make sure to set the "display" option to off, otherwise you get 11^4 notifications.

This took about 60 seconds.

The one thing that I did not check was wether the way you have the anonymous function set up correctly.

I wonder if Lx(i) and Ly(j)... etc. etc. in the function is changing.

There are several ways to check but I would double check that they are...

Although this works, that function is painfully ugly... There has got to be a way to clean this up.

clc
clear
tic
fy=40;
cover=1.25;
ql=1;
Lx = zeros(11,1);
Ly = zeros(11,1);
Cx = zeros(11,1);
Cy = zeros(11,1);
B = zeros(11,11);
A = zeros(11,11,11,11);
D = zeros(11,11,11,11);
oldD = 3.;
options = optimset('Display','off');
for k=1:11
%k = 1;
    Lx(k)=120+(k-1)*24;
for l=1:11
%l = 1;
    Ly(l)=120+(l-1)*24;
for i=1:11
%i = 1;
    Cx(i)=8+(i-1)*4;
for j=1:11
%j = 1;
    Cy(j)=8+(j-1)*4;

if (Cx(i)>=Cy(j))
    B(i,j)=Cx(i)/Cy(j);
else B(i,j)=Cy(j)/Cx(i);
end
A(i,j,k,l)=(Lx(k)+Ly(l))/(Cx(i)+Cy(j));

newD = fsolve(@(d)((1+((5/32)*((d*(Cx(i)+Cy(j)+2*d)/((Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d)) ...
    *(1+(3/2)*sqrt((Cy(j)+d)/(Cx(i)+d))))*0.336*Ly(l)*(Lx(k)-Cx(i))^2/ ...
    (d^3+(Cx(i)+d)^2*d+3*d*(Cy(j)+d)*(Cx(i)+d)))+(d*(Cx(i)+Cy(j)+2*d)/ ...
    ((Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))*(1+(3/2)*sqrt((Cx(i)+d)/(Cy(j)+d)))) ...
    *0.336*Lx(k)*(Ly(l)-Cy(j))^2/(d^3+(Cy(j)+d)^2*d+3*d*(Cx(i)+d)*(Cy(j)+d)))))) ...
    *6.4*ql*(Lx(k)*Ly(l)-(Cx(i)+d)*(Cy(j)+d))/(d*(2*Cx(i)+2*Cy(j)+4*d)) ...
    -min([2+4/B(i,j);(40*d)/((2*(Cx(i)+d)+2*(Cy(j)+d)))+2;4])), oldD, options);
D(i,j,k,l) = newD;
oldD = newD;
end
end
end
end

toc