Some advice after taking a look at this code.
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Hornett
le 30 Sep 2024
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
0 commentaires
Voir également
Catégories
En savoir plus sur Data Analysis dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!