Some advice after taking a look at this code.

Question

Ricardo Boza Villar 2022-4-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1701045-some-advice-after-taking-a-look-at-this-code

回答： Nihal 2024-9-30，11:23

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Nihal 2024-9-30，11:23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1701045-some-advice-after-taking-a-look-at-this-code#answer_1524710

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