Some advice after taking a look at this code.

Hello, I want to improve. I welcome suggestions. Thanks.

Main way to run is with VigaOrden4(1,@(x) 2 + sin(2*pi*x)) although the solution is rather anticlimactic because the deformation is so small.

function VigaOrden4(kappa, funr)
%% Ejemplos
% VigaOrden4(1,@(x) 2 + sin(2*pi*x))
% VigaOrden4(1,@(x) cos(x))
% VigaOrden4(1,@(x)1 + sin(x))
%% Datos
k = kappa;
r = funr;
h = 1e-5;
%% Condiciones de contorno
bcfun = @(ua,ub)[ua(1); ua(2); ub(1); ub(2)]; % no estoy seguro de cuáles son las condiciones de contorno ni de cómo escribirlas
%% Inicializaciones
xinit = linspace(0,1,10);
yinit = [0; 0; 0; 0];
solinit = bvpinit(xinit, yinit);
%% Opciones bvp4c
options = bvpset('RelTol',1e-6);
%% Resolución numérica con bvp4c
sol = bvp4c(@(x, u)f(x, u, k, r, h), bcfun, solinit, options);
% length(sol.x)
%% Interpolación de la solución
x = linspace(0,1,1e3);
u = deval(sol, x, 1);
%% Dibujo de la solución
close all
%% Figura 1
figure(1), hold on
plot(x, u, 'b', 'LineWidth', 2)
axis([0 1 -3e-3 1e-3])
title('Vector deformación')
%% Figura 2
figure(2), hold on
plot(x, u, 'b', 'LineWidth', 2)
axis([0 1 -3e-3 1e-3])
axis equal
title('Vector deformación en perspectiva')
%% Figura 3
figure(3), hold on, dibuja(r, zeros(1,length(u)), x), view([30 25]), title('Sin deformación')
%% Figura 4
figure(4), hold on, dibuja(r, u, x), view([30 25]), title('Con deformación')
%% Figura 5 (porque es la última orden, puedo machacar variables)
epsilon = 1;
x = linspace(0,1);
y = linspace(-2,2);
[vx, vy] = meshgrid(x,y);
z = real(sqrt(epsilon*r(vx) - vy.^2));
figure(5), hold on
surf(vx,vy,z)
surf(vx,vy,-z)
view([30 25])
axis equal
%% Función principal del problema de contorno
    function [dudt] = f(x, u, k, r, h)
        dudt = [u(2); u(3); u(4); -1/(k*r(x)^2)-8*(rp(r,x,h)/r(x))*u(4)-4*(3*(rp(r,x,h)/r(x))^2+rpp(r,x,h)/r(x))*u(3)];
    end
%% Derivada primera función r (aproximación numérica)
    function [drdt] = rp(r, x, h)
        drdt = (r(x+h) - r(x))/h;
    end
%% Derivada segunda función r (aproximación numérica)
    function [d2rdt2] = rpp(r, x, h)
        d2rdt2 = (r(x+h) + r(x-h) -2*r(x))/h^2;
    end
%% Dibujo de la superficie del sólido antes de la deformación
    function dibuja(r,u,x)
        n = 50;
        xx = linspace(0,1,n);
        y = linspace(-2,2,n);
        newu = interp1(x,u,xx); % error muy burdo, pero bueno. Todo sea por hacer la gráfica
        % whos newu
        for m = 1:5:n
            z = real(sqrt(r(xx(m))-y.^2)-newu(m)); % -u
            xxx = xx(m)*ones(n);
            plot3(xxx, y, z,'b')
            plot3(xxx, y, -z,'b')
        end
    end
end

Answer 1

[Hornett]

Hi Ricardo,

After going through your code i would suggest you to avoid unnecessary function calls, Minimize the number of function calls within loops. For instance, you can calculate rp(r, x, h) and rpp(r, x, h) once and store the results in variables before using them.

I hope it helps