Can anyone help me to fix my code?

Can someone please tell me what am I doing wrong in following code please?

clear all
clc
syms x y a b 
 
X=[1:1:10];
Y=[0.5379,0.827,1.2654,1.5324,1.717,1.3983,1.6844,1.6892,1.506,1.707];
 
 
a = [1;0;0;0;0;0;0;0;0;0;0];
b = [1;0;0;0;0;0;0;0;0;0;0];
c = [0;0;0;0;0;0;0;0;0;0;0];
 
u = [((a(1)*x)/(b(1)+x)) , ((a(1)*2*x)/(b(1)+2*x)), ((a(1)*3*x)/(b(1)+3*x)), ((a(1)*4*x)/(b(1)+4*x)), ((a(1)*5*x)/(b(1)+5*x)), ((a(1)*6*x)/(b(1)+6*x)), ((a(1)*7*x)/(b(1)+7*x)), ((a(1)*8*x)/(b(1)+8*x)), ((a(1)*9*x)/(b(1)+9*x)), ((a(1)*10*x)/(b(1)+10*x))]';
 
B0=[a(1);b(1)];
U = [0;0;0;0;0;0;0;0;0;0;0];
 
 
for i=2:11
    
    U(i-1) = [((a(i-1)*1)/(b(i-1)+1)) , ((a(i-1)*2)/(b(i-1)+2)), ((a(i-1)*3)/(b(i-1)+3)), ((a(i-1)*4)/(b(i-1)+4)), ((a(i-1)*5)/(b(i-1)+5)), ((a(i-1)*6)/(b(i-1)+6)), ((a(i-1)*7)/(b(i-1)+7*x)), ((a(i-1)*8)/(b(i-1)+8)), ((a(i-1)*9)/(b(i-1)+9)), ((a(i-1)*10)/(b(i-1)+10))]';
 
    Jac = jacobian([((a*1)/(b+1)) , ((a*2)/(b+2)), ((a*3)/(b+3)), ((a*4)/(b+4)), ((a*5)/(b+5)), ((a*6)/(b+6)), ((a*7)/(b+7)), ((a*8)/(b+8)), ((a*9)/(b+9)), ((a*10)/(b+10))],  [a,b] );
 
    r = u - U(i-1);
    
    c = [a(i-1);b(i-1)] - inv(J(i-1)'.* J(i-1)).*(J(i-1)'.*r(i-1));
    
    c
    
    a(i) = c(1);
    b(i) = c(2);
end

Answer 1

[David Goodmanson]

Hello JCam,

Since you are looking for a numerical result there is really no reason to use syms here. With x and y both being 10x1 column vectors. the basic idea is to minimize [y - a*x/((b+x)]. SInce there are 10 equations and two unknowns, this is going to be done in the least squares sense.

The basic approach for some function f(z) that is suppossed to equal a target value is

f(z) + df/dz(z)*(znew - z) = target value  or
       df/dz(z)*(znew - z) = target value - f(z)     

For variables a and b the Jacobian of the fitting function is

d/da [a*x/(b+x)] = [x/(b+x)]
d/db [a*x/(b+x)] = [-a*x/(b+x)^2]

where everything in brackets represents 10x1 column vectors. Let J be the 10x2 matrix constructed by putting the two known column vectors above side-by-side. Let B = [a;b] be a 2x1 matrix. Then for the iteration in matrix notation,

J*(B_n+1 - B_n) = [y - a*x/((b+x)]   where a = B_n(1) and b = B_n(2).

This is what the notes have, except for some reason they define the Jacobian with a minus sign and then they use a compensating minus sign in what they call the recursive equation. (It's really an iterated equation, not recursive). The code below does not have the minus signs.

To get (B_n+1 - B_n) you have to do something about J. The approach in the notes is to multiply both sides of the equation above by J', producing J'*J which is 2x2. Then you multiply both sides by the inverse of that matrix for the final result. That's the so-called normal equation approach, which is great if it were still 1950. Today you can use Matlab left divide, which is faster and more accurate:

(B_n+1 - B_n) = J\[y - a*x/((b+x)]

(since J'*J is 2x2 it probably does not matter numerically to use the normal equation here, but there is no reason to promote an obsoete method. It's fun to go riding in a 1950 Buick convertible, but inferior dated algorithms don't have the same appeal).

% column vectors
x = (1:1:10)';
y = [0.5379,0.827,1.2654,1.5324,1.717,1.3983,1.6844,1.6892,1.506,1.707]';
B = [1 1]'
for k = 1:10
 a = B(1);
 b = B(2);
 J  = [x./(b+x) -a*x./(b+x).^2];
 Bnew = B + J\(y - a*x./(b+x))
 B = Bnew;
end  
plot(x,y,x,a*x./(b+x))
grid on