Cannot get a non-singular matrix

Hello all, I'm trying to make a non-singular matrix, but I don't know how to solve this problem.

If I use 11 instead of 13 as the value of m, I get a non-singular matrix.

f = @(x) log(exp(x) + 1) + 2 * sin(x);
a1 = 1;
b1 = 18;
m = 13;
x = linspace(a1, b1, m)';
y = f(x);
b = @(x, j) x.^(j - 1);
B = ones(13,13);
for j = 2 : m
    B(:, j) = b(x, j);
end

a = B\y;

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.034785e-19.

Answer 1

[Sajid Afaque]

you need to focus on the last step i.e. division part

a = B ./ y ; %element wise division

[Giuseppe Fiorillo]

But using the ./ operator I will get the result of the division between B and y, but I need to get the solution to the sistem Ba = y, so I need the backslash (\) operator.

Answer 2

[Mitch Lautigar]

Couple comments. Please see my added comments to the code

f = @(x) log(exp(x) + 1) + 2 * sin(x); %I'd change "*" to ".*" since you are passing in an array

a1 = 1;

b1 = 18;

m = 13;

x = linspace(a1, b1, m)'; %this line of code works fine. I tested it

y = f(x); %if you make the change suggested above, should run through clearly.

b = @(x, j) x.^(j - 1); %I'm hoping j is initialized elsewhere.

B = ones(13,13); %I'd make this line read "B = ones(m,m);" unless you need a static . %13x13 array

for j = 2 : m

B(:, j) = b(x, j);

end

a = B\y;

[Giuseppe Fiorillo]

I tried changing "*" to ".*" but the matrix remains singular.

j isn't initialized elsewhere because it's a parameter of b (I will use it anyway in the for cicle).

B=ones(m, m) is correct!

Answer 3

[Torsten]

Maybe of help: The Vandermonde matrix is ill conditioned !

Maybe of help:

f = @(x) log(exp(x) + 1) + 2 * sin(x);
a1 = 1;
b1 = 18;
m = 13;
x = linspace(a1, b1, m)';
y = f(x);
b = @(x, j) x.^(j - 1);
B = ones(m,m);
for j = 2 : m
    B(:, j) = b(x, j);
end
a = B\y;

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.034785e-19.

xx = linspace(a1,b1,100).';
yapprox = zeros(numel(xx),1);
for i=1:numel(xx)
    summe = 0.0;
    for j=1:m
        summe = summe + a(j)*xx(i)^(j-1);
    end
    yapprox(i) = summe;
end
plot(xx,f(xx))
hold on
plot(xx,yapprox)