sqp() with trapz() Classic List Threaded 17 messages Open this post in threaded view
|

sqp() with trapz()

 Hello all, I have been working for a while in a code in Octave. All lines seem to work, except for the last one, when I run an optimization using sqp(). This is the code: %Use Gnuplot for graphics in Octave graphics_toolkit("gnuplot") %Exact solution F=@(mu) mu*log((-1-sqrt(1+mu^2))/(1-sqrt(1+mu^2)))-1; %Function for the solver mu=fzero(F,[0.01,1]) %Solution via Brent-Dekker-van Wijngaarden method s=linspace(0,2,N=51); %s discretization y=-sqrt(1+mu^2)+sqrt((s-1).^2+mu^2); %y exact values %plot(s,y) %Exact solution %Approximate solution n=50; %Number of subintervals h=2/n; %Step size F=@(y) trapz(s,y); %Approximate objective function (Trapezoidal rule) D=zeros(n+1,n+1); D1=-eye(n/2)+diag(ones([n/2-1,1]),1); D4=eye(n/2)+diag(-ones([n/2-1,1]),-1); D(1:n/2,1:n/2)=D1; D(end-n/2+1:end,end-n/2+1:end)=D4; D(n/2,n/2+1)=1; D(n/2+1,n/2+2)=1/2; D(n/2+2,n/2+1)=-1; D(n/2+1,n/2)=-1/2; D=1/h*D; %Differentiation matrix yp=@(y) transpose(D*transpose(y)); %Derivative (Finite Differences) G=@(y) sqrt(1-yp(y).^2); %Function in the integral constraint R=@(y) trapz(s,G(y))-1; %Approximation of the integral constraint (Trapezoidal rule) sqp(y,F,R) %Optimization via sequential quadratic programming And the error it returns: error: trapz: X and Y must have same shape error: called from     trapz at line 120 column 9     sqp at line 351 column 7     Catenary at line 28 column 1 However the trapz() function alone seems to work correctly. Does anyone know what is happening here? Thank you in advance for any help! -- Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
Open this post in threaded view
|

Re: sqp() with trapz()

 Am 25.08.19 um 17:18 schrieb dg.aragones: > Hello all, > > I have been working for a while in a code in Octave. All lines seem to work, > except for the last one, when I run an optimization using sqp(). This is the > code: > > %Use Gnuplot for graphics in Octave > graphics_toolkit("gnuplot") > > %Exact solution > F=@(mu) mu*log((-1-sqrt(1+mu^2))/(1-sqrt(1+mu^2)))-1; %Function for the > solver > mu=fzero(F,[0.01,1]) %Solution via Brent-Dekker-van Wijngaarden method > s=linspace(0,2,N=51); %s discretization > y=-sqrt(1+mu^2)+sqrt((s-1).^2+mu^2); %y exact values > %plot(s,y) %Exact solution > > %Approximate solution > n=50; %Number of subintervals > h=2/n; %Step size > F=@(y) trapz(s,y); %Approximate objective function (Trapezoidal rule) > D=zeros(n+1,n+1); > D1=-eye(n/2)+diag(ones([n/2-1,1]),1); > D4=eye(n/2)+diag(-ones([n/2-1,1]),-1); > D(1:n/2,1:n/2)=D1; > D(end-n/2+1:end,end-n/2+1:end)=D4; > D(n/2,n/2+1)=1; > D(n/2+1,n/2+2)=1/2; > D(n/2+2,n/2+1)=-1; > D(n/2+1,n/2)=-1/2; > D=1/h*D; %Differentiation matrix > yp=@(y) transpose(D*transpose(y)); %Derivative (Finite Differences) > G=@(y) sqrt(1-yp(y).^2); %Function in the integral constraint > R=@(y) trapz(s,G(y))-1; %Approximation of the integral constraint > (Trapezoidal rule) > sqp(y,F,R) %Optimization via sequential quadratic programming > > And the error it returns: > > error: trapz: X and Y must have same shape > error: called from >      trapz at line 120 column 9 >      sqp at line 351 column 7 >      Catenary at line 28 column 1 > > However the trapz() function alone seems to work correctly. Does anyone know > what is happening here? > > Thank you in advance for any help! > > > > -- > Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html> > Did you look up the description of trapz? https://octave.sourceforge.io/octave/function/trapz.html
Open this post in threaded view
|

Re: sqp() with trapz()

 In reply to this post by dg.aragones On Sun, Aug 25, 2019 at 11:18 AM dg.aragones <[hidden email]> wrote:Hello all, I have been working for a while in a code in Octave. All lines seem to work, except for the last one, when I run an optimization using sqp(). This is the code: %Use Gnuplot for graphics in Octave graphics_toolkit("gnuplot") %Exact solution F=@(mu) mu*log((-1-sqrt(1+mu^2))/(1-sqrt(1+mu^2)))-1; %Function for the solver mu=fzero(F,[0.01,1]) %Solution via Brent-Dekker-van Wijngaarden method s=linspace(0,2,N=51); %s discretization y=-sqrt(1+mu^2)+sqrt((s-1).^2+mu^2); %y exact values %plot(s,y) %Exact solution %Approximate solution n=50; %Number of subintervals h=2/n; %Step size F=@(y) trapz(s,y); %Approximate objective function (Trapezoidal rule) D=zeros(n+1,n+1); D1=-eye(n/2)+diag(ones([n/2-1,1]),1); D4=eye(n/2)+diag(-ones([n/2-1,1]),-1); D(1:n/2,1:n/2)=D1; D(end-n/2+1:end,end-n/2+1:end)=D4; D(n/2,n/2+1)=1; D(n/2+1,n/2+2)=1/2; D(n/2+2,n/2+1)=-1; D(n/2+1,n/2)=-1/2; D=1/h*D; %Differentiation matrix yp=@(y) transpose(D*transpose(y)); %Derivative (Finite Differences) G=@(y) sqrt(1-yp(y).^2); %Function in the integral constraint R=@(y) trapz(s,G(y))-1; %Approximation of the integral constraint (Trapezoidal rule) sqp(y,F,R) %Optimization via sequential quadratic programming And the error it returns: error: trapz: X and Y must have same shape error: called from     trapz at line 120 column 9     sqp at line 351 column 7     Catenary at line 28 column 1 However the trapz() function alone seems to work correctly. Does anyone know what is happening here? Thank you in advance for any help! -- Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html You have 2 different functions named Fwhich one do you want in the aqp function?-- DAS Open this post in threaded view
|

Re: sqp() with trapz()

 In reply to this post by karl Thank you for your answer! Yes, I did. However, I think I've written the code accordingly... -- Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
Open this post in threaded view
|

Re: sqp() with trapz()

 In reply to this post by Doug Stewart-4 Thank you for your answer! It's true. I want to use the second one in the sqp() function. -- Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
Open this post in threaded view
|

Re: sqp() with trapz()

 In reply to this post by dg.aragones Am 26.08.19 um 01:06 schrieb dg.aragones: > Thank you for your answer! > > Yes, I did. However, I think I've written the code accordingly... > > > > -- > Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html> > The error message says that         error: trapz: X and Y must have same shape You have written: R=@(y) trapz(s,G(y))-1 Do s and G(y) have the same shape?
Open this post in threaded view
|

Re: sqp() with trapz()

 Yes, they have the same size. It's weird... -- Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
Open this post in threaded view
|

Re: sqp() with trapz()

 Am 26.08.19 um 18:25 schrieb dg.aragones: > Yes, they have the same size. It's weird... > > > > -- > Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html> > Please look at your code. Trapz expects two vectors of the same size. You have the placeholder variable s (replace it by x since it can any name) and then a function call. Are these two vectors of the same size?
Open this post in threaded view
|

Re: sqp() with trapz()

 Am 26.08.19 um 19:07 schrieb [hidden email]: > Am 26.08.19 um 18:25 schrieb dg.aragones: >> Yes, they have the same size. It's weird... >> >> >> >> -- >> Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html>> >> > Please look at your code. Trapz expects two vectors of the same size. > > You have the placeholder variable s (replace it by x since it can any name) > > and then a function call. Are these two vectors of the same size? > > > Replace s by y in the answer.
Open this post in threaded view
|

Re: sqp() with trapz()

 Yes, they have exactly the same size: s is a vector 51x1 G(y) returns a vector 51x1 -- Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
Open this post in threaded view
|

Re: sqp() with trapz()

 Am 26.08.19 um 20:12 schrieb dg.aragones: > Yes, they have exactly the same size: > > s is a vector 51x1 > G(y) returns a vector 51x1 > > > > -- > Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html> > if you mean, I give up.
Open this post in threaded view
|

Re: sqp() with trapz()

 In reply to this post by dg.aragones On 8/26/19 2:12 PM, dg.aragones wrote: > Yes, they have exactly the same size: > > s is a vector 51x1 > G(y) returns a vector 51x1 > The error happens in the other invocation of trapz, via the F function F=@(y) trapz(s,y);
Open this post in threaded view
|

Re: sqp() with trapz()

 Yes, but I don't understand why, because they have also the same size. They are both 51x1 vectors. -- Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html
Open this post in threaded view
|

Re: sqp() with trapz()

 On 8/26/19 2:32 PM, dg.aragones wrote: > Yes, but I don't understand why, because they have also the same size. They > are both 51x1 vectors. Nope, one is 1x51 and the other is 51x1. By using @ functions, it's harder to debug this because you can't put tests/printouts into them, so I redefined them to use the function trapz_debug, and defined it like so: function f = trapz_debug (a, b)    size (a)    size (b)    f = trapz (a, b) endfunction
Open this post in threaded view
|

Re: sqp() with trapz()

 Ok, thank you very much for your answer! -- Sent from: https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html