Some advice after taking a look at this code.
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/971015/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/971020/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/971025/image.png)
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 Kommentare
Antworten (1)
Nihal
am 10 Nov. 2023
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 Kommentare
Siehe auch
Kategorien
Mehr zu Downloads finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!