Symbolic Substitution in While Loop - What am I missing?
in [Matlab]
|
From: Crystal Nassouri on 30 Mar 2010 14:18 Hi there, I am having an issue with my while loop. I have a symbolic equation that gets values further down the program, but after the first iteration of the while loop, it seems to reset to the symbolic values instead of holding the numerical values. My code is below. I am used to C style programming and am having trouble figuring out how to code similarly in Matlab. Thanks in advance! Crystal %Prob 4.2.2 clear all; syms x1; syms x2; syms p; syms Lk; syms c; n = 0; optL = -2; Lnext = 0; lamda0 = 0; p = 1; ck = 1; %for all k, c = 1 %Lc = (1/2)*(x1^2+ p*(x2^p))+ 2*x2 + lamda0*x2+(ck/2)*(x2^2); Lc = (1/2)*(x1^2-x2^2)-3*x2+Lk*x2+(c/2)*x2^2; gx1 = diff(Lc,x1) gx2 = diff(Lc,x2) x1k = solve(gx1,x1) x2k = solve(gx2,x2) %Unable to solve for x2 explicity. Also having %difficulty with some p's H = x2k Lk = lamda0 while Lnext > -2 c = 10^(n+1) Lnext = Lk + c*H %ARRRG WHY ISNT IT SUBSTITUTING? Lk = Lnext n = n+1; end disp(['Iterations for convergence: n=' num2str(n)]);
From: Walter Roberson on 30 Mar 2010 14:53 Crystal Nassouri wrote: Lc = (1/2)*(x1^2-x2^2)-3*x2+Lk*x2+(c/2)*x2^2; > gx1 = diff(Lc,x1) > gx2 = diff(Lc,x2) > x1k = solve(gx1,x1) > x2k = solve(gx2,x2) %Unable to solve for x2 explicity. Also having > %difficulty with some p's Huh? Solving gx2 for x2 is trivial. Distribute the 1/2 over the sum in Lc to get 1/2*x1^2 - 1/2*x2^2 - 3*x2 + Lk * x2 + (c/2)*x2^2 . Now differentiate with respect to x2: 0 - 2*(1/2)*x2 - 3 + Lk + 2*(c/2)*x2 which is -x2 - 3 + Lk + c*x2 which is (c-1)*x2 -3 + Lk Solve for 0 trivially: x2 = (3 - Lk)/(c-1) However, if you use your original Lc value, commented out, of Lc = (1/2)*(x1^2+ p*(x2^p))+ 2*x2 + lamda0*x2+(ck/2)*(x2^2); then Yes, you are going to have problems. If p is 6 or higher, or p is not a positive integer, then after the differentiation you are going to be left with an expression which is not a polynomial of degree 4 or less, and by the Fundamental Theorem of Algebra, you cannot analytically find the roots of anything higher than a quartic polynomial. > H = x2k > Lk = lamda0 > > while Lnext > -2 > c = 10^(n+1) > Lnext = Lk + c*H %ARRRG WHY ISNT IT SUBSTITUTING? > Lk = Lnext > n = n+1; > end If you want to force substitution, then use subs() around the symbol or expression.
From: Crystal Nassouri on 30 Mar 2010 16:05 Hi, Thanks for your reply. Even though solving gx2 for x2 might be trivial, I am trying to write out a general case program for iteratively solving for an optimal Lagrangian. The problem is in: > while Lnext > -2 > c = 10^(n+1) > Lnext = Lk + c*H %ARRRG WHY ISNT IT SUBSTITUTING? > Lk = Lnext > n = n+1; > end Where Lnext does take in the val for Lk and c, but does not take H numerically. When I run the program, I can see that Lnext is still symbolic because for some reason H is symbolic even though I have given it all required values. This kind of thing would run in C or C++, do I have to "force" substitution every time I want to do something like this? I'd like to start using Matlab for things like this but at this point I am spending more time trying to get it to do simple things when I could have written this and more in C. "Crystal Nassouri" <crystal.nassouri(a)gmail.com> wrote in message <hotf9d$br7$1(a)fred.mathworks.com>... > Hi there, > I am having an issue with my while loop. I have a symbolic equation that gets values further down the program, but after the first iteration of the while loop, it seems to reset to the symbolic values instead of holding the numerical values. My code is below. I am used to C style programming and am having trouble figuring out how to code similarly in Matlab. > > Thanks in advance! > Crystal > > %Prob 4.2.2 > > clear all; > syms x1; > syms x2; > syms p; > syms Lk; > syms c; > n = 0; > optL = -2; > Lnext = 0; > > lamda0 = 0; > p = 1; > > ck = 1; %for all k, c = 1 > > %Lc = (1/2)*(x1^2+ p*(x2^p))+ 2*x2 + lamda0*x2+(ck/2)*(x2^2); > Lc = (1/2)*(x1^2-x2^2)-3*x2+Lk*x2+(c/2)*x2^2; > > gx1 = diff(Lc,x1) > gx2 = diff(Lc,x2) > > x1k = solve(gx1,x1) > x2k = solve(gx2,x2) %Unable to solve for x2 explicity. Also having > %difficulty with some p's > H = x2k > Lk = lamda0 > > while Lnext > -2 > c = 10^(n+1) > Lnext = Lk + c*H %ARRRG WHY ISNT IT SUBSTITUTING? > Lk = Lnext > n = n+1; > end > > disp(['Iterations for convergence: n=' num2str(n)]); > > >
From: Walter Roberson on 30 Mar 2010 16:20 Crystal Nassouri wrote: > The problem is in: >> while Lnext > -2 >> c = 10^(n+1) >> Lnext = Lk + c*H %ARRRG WHY ISNT IT SUBSTITUTING? >> Lk = Lnext >> n = n+1; >> end > > Where Lnext does take in the val for Lk and c, but does not take H > numerically. When I run the program, I can see that Lnext is still > symbolic because for some reason H is symbolic even though I have given > it all required values. If you have a numeric symbolic expression that you want to evaluate to a hardware floating point number, call double() on the expression and save the result. In your case you could do that outside the loop (because H does not change within the loop.)
From: Steven Lord on 30 Mar 2010 17:57
"Walter Roberson" <roberson(a)hushmail.com> wrote in message news:hotmeu$kji$2(a)canopus.cc.umanitoba.ca... > Crystal Nassouri wrote: > >> The problem is in: >>> while Lnext > -2 >>> c = 10^(n+1) >>> Lnext = Lk + c*H %ARRRG WHY ISNT IT SUBSTITUTING? >>> Lk = Lnext >>> n = n+1; >>> end >> >> Where Lnext does take in the val for Lk and c, but does not take H >> numerically. When I run the program, I can see that Lnext is still >> symbolic because for some reason H is symbolic even though I have given >> it all required values. > > If you have a numeric symbolic expression that you want to evaluate to a > hardware floating point number, call double() on the expression and save > the result. In your case you could do that outside the loop (because H > does not change within the loop.) You can't substitute into H outside the loop, as H depends on both Lk and c which _are_ changing inside the loop. The OP can use SUBS to substitute values into H or convert H into a function handle of c and Lk using matlabFunction and then simply evaluate it as hHandle(c, Lk) inside the loop. If you do the latter, though, beware when n gets too large -- you could lose precision because 10^q cannot be exactly represented in double when q > 15. -- Steve Lord slord(a)mathworks.com comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ |