convergent with order p = 1

I am trying to write a program that produces a figure of convergence with order 𝑝 = 1.

But I get an error message on polyfit function that says:

Error using polyfit

The first two inputs must have the same number of elements.

Would you please help me to figure what is the problem?

The function is as following:

% Creating a function of semi-implicit Euler method
function u = semi_imp_euler(x0, v0, nsteps, tend, k, m)
dt = tend/double(nsteps);
u = zeros(2, nsteps+1);
u(1,1) = x0;
u(2,1) = v0;
for n = 1:nsteps
    u(:,n+1) = [(1-(dt^(2))*(k/m))*u(1,n)+dt*u(2,n); (-dt*(k/m))*u(1,n)+u(2,n)];
end
end

---------------------

The program that I am getting a message from:

T = 10.0; % Final time until which we compute
N = [1000, 750, 500, 250, 100, 75, 50, 10]; % Number of time steps
k = 5;
m = 0.5;
% Initial values
x0 = 1.0;
v0 = 0.1;
alpha = sqrt(k/m); % Computing alpha for linear solution
% Computing parameters A and B for exact linear solution 
A = v0/alpha;
B = x0;
% Computing exact linear solution
u_exact = @(t) [A*sin(alpha*t)+B*cos(alpha*t); v0*cos(alpha*t)-x0*alpha*sin(alpha*t)];
error_semi_imp = zeros(2,length(N));
dts = zeros(1,length(N));
for n = 1:length(N)
    u_semi_imp_euler = semi_imp_euler(x0, v0, N(n), T, k, m);
    taxis = linspace(0, T, N(n)+1);
    dts(n) = taxis(2) - taxis(1);  % store the time step dt for plotting
    error_semi_imp(:,n) = max(abs(u_semi_imp_euler - u_exact(taxis)),[],2);
end
p_semi_imp = polyfit(log(dts), log(error_semi_imp), 1);

---------------

The error message:

Error using polyfit

The first two inputs must have the same number of elements.

Error in

p_semi_imp = polyfit(log(dts), log(error_semi_imp), 1);

Answer 1

[Walter Roberson]

T = 10.0; % Final time until which we compute
N = [1000, 750, 500, 250, 100, 75, 50, 10]; % Number of time steps
k = 5;
m = 0.5;
% Initial values
x0 = 1.0;
v0 = 0.1;
alpha = sqrt(k/m); % Computing alpha for linear solution
% Computing parameters A and B for exact linear solution 
A = v0/alpha;
B = x0;
% Computing exact linear solution
u_exact = @(t) [A*sin(alpha*t)+B*cos(alpha*t); v0*cos(alpha*t)-x0*alpha*sin(alpha*t)];
error_semi_imp = zeros(2,length(N));
dts = zeros(1,length(N));
for n = 1:length(N)
    u_semi_imp_euler = semi_imp_euler(x0, v0, N(n), T, k, m);
    taxis = linspace(0, T, N(n)+1);
    dts(n) = taxis(2) - taxis(1);  % store the time step dt for plotting
    error_semi_imp(:,n) = max(abs(u_semi_imp_euler - u_exact(taxis)),[],2);
end

size(dts)

ans = 1×2

1 8

size(error_semi_imp)

ans = 1×2

2 8

p_semi_imp = polyfit(log(dts), log(error_semi_imp), 1);

Error using polyfit
The first two inputs must have the same number of elements.

% Creating a function of semi-implicit Euler method
function u = semi_imp_euler(x0, v0, nsteps, tend, k, m)
    dt = tend/double(nsteps);
    u = zeros(2, nsteps+1);
    u(1,1) = x0;
    u(2,1) = v0;
    for n = 1:nsteps
        u(:,n+1) = [(1-(dt^(2))*(k/m))*u(1,n)+dt*u(2,n); (-dt*(k/m))*u(1,n)+u(2,n)];
    end
end

Look at your semi_imp_euler code: it is hard-coded to initialize u(2,1) so the output will have at least two rows -- and each step you set u(:,n+1) to a vector of 2 elements, so u will certainly have two rows on output.

You then take both rows into u_semi_imp_euler and you max() a calculation with it along dimension 2, so the result of the max() is going to be 2 x 1. So error_semi_imp is going to be constructed as 2 rows and n columns.

Then you try to polyfit something of length n (dts) as the independent variable, and something of size 2 x n as the dependant variable. But the independent variable for polyfit must be a vector the same length as the dependent variable.

If you want to polyfit the two rows of error_semi_imp independently then you need to make two polyfit calls.

[Abdullah]

How can I make two polyfit calls for the two rows?

[Walter Roberson]

p_semi_imp_x = polyfit(log(dts), log(error_semi_imp(1,:)), 1);
p_semi_imp_v = polyfit(log(dts), log(error_semi_imp(2,:)), 1);

[Abdullah]

Thanks. This worked for the polyfit, but the result for each is not for order p = 1.

The program produced the figures as following:

I can see the points on vertical line 10^(0) in both figures that might cause the non-linear answer, but I do not know why this happens.

My code for plotting:

p_semi_imp_x = polyfit(log(dts), log(error_semi_imp(1,:)), 1);
p_semi_imp_v = polyfit(log(dts), log(error_semi_imp(2,:)), 1);

figure(1);
loglog(dts, error_semi_imp(1,:), 'ro'); hold on;
loglog(dts, exp(polyval(p_semi_imp_x, log(dts))), 'r-');
txt = strcat('Slope p=', num2str(p_semi_imp_x(1),'%3.2f'));
text(1e-2, 1e-3, txt); 
xlim([dts(1) dts(end)]);
xlabel('\Delta t');
ylabel('Error');
grid on;
legend('semi_imp_x','Linear Fit_x')

figure(2);
loglog(dts, error_semi_imp(2,:), 'bo'); hold on;
loglog(dts, exp(polyval(p_semi_imp_v, log(dts))), 'b-');
txt = strcat('Slope p=', num2str(p_semi_imp_v(1),'%3.2f'));
text(1e-2, 1e-3, txt); 
xlim([dts(1) dts(end)]);
xlabel('\Delta t');
ylabel('Error');
grid on;
legend('semi_imp_v','Linear Fit_v')

Could you please help me to figure out the problem?