Bonsoir je voudrai savoir comment tracé ces fonctions du code ci-joint sur matlab ? Merci beaucoup
voila le code: Code: function [L2_error_w] = beam_loadsteps(N,LS) % problem parameters param.L = 1; param.EA = 10^5; param.EI = 20; param.p1 = 0; param.p3 = 10; param.num_load_steps = linspace(0,1,LS); F1crit=403.814571; F1=0.50*F1crit; % problem kinematic and static boundary conditions bc.u_idx = [1 ]; bc.u_val = [0 ]; bc.w_idx = [2 3*N-1]; bc.w_val = [0 0]; bc.phi_idx = [3]; bc.phi_val = [0]; bc.F1_idx = [3*N-2]; bc.F1_val = [-F1]; bc.F3_idx = [ ]; bc.F3_val = [ ]; bc.M_idx = [ ]; bc.M_val = [ ]; % define number of points for evaluating exact solution E = N*100; % obtain exact and approximate solution %[x_exact ,U_exact ] = beam_exact(E, param, bc); [x_exact ,U_exact ] = beam_exact_loadsteps(E, param, bc); [x_approx,U_approx] = beam_approx_loadsteps(N, param, bc); % prepare plotting close all hidden; figure('DefaultAxesFontSize',18); clf; hold on; set(gca,'Ydir','reverse'); % get sub-element values (for plotting and error computation) u_exact = U_exact(1,:); w_exact = U_exact(2,:); u_approx_m = interp1 (x_approx,U_approx(1,:,end), x_exact); w_approx_m = interp1hermite(x_approx,U_approx(2,:,end),U_approx(3,:,end),x_exact); % plot exact and approximate solution graph plot(x_exact ,u_exact , 'b--'); plot(x_exact ,u_approx_m, 'b-' ); plot(x_exact ,w_exact , 'r--'); plot(x_exact ,w_approx_m, 'r-' ); plot(x_approx,U_approx(1,:,end),'bo'); plot(x_approx,U_approx(2,:,end),'ro'); axis([0 param.L -0.005 +0.015]); xlabel('coordinate x [m]'); ylabel('displacement [m]'); title('beam 2nd order theory'); legend('u_{exact}^{2nd order}', 'u_{approx}^{2nd order}', 'w_{exact}^{2nd order}', 'w_{approx}^{2nd order}', 'Location','SouthWest'); grid on; print('beam_2ndorder.png','-dpng'); % plot lambda-deformation curve figure('DefaultAxesFontSize',18); clf; hold on; phi = squeeze(U_approx(3,end,:)); plot(phi, param.num_load_steps, 'b-*'); xlabel('deformation value (end point rotation) [rad]'); ylabel('load factor [-]'); set(gca,'Xdir','reverse'); title('beam 2nd order theory: load-deformation curve'); grid on; print('beam_2ndorder_loadcurve.png','-dpng'); % compute L_2 error w L2_error_w = sqrt( trapz( (w_approx_m - w_exact).^2 ) / trapz(w_exact.^2) ); end
function [L2_error_w] = beam_loadsteps(N,LS) % problem parameters param.L = 1; param.EA = 10^5; param.EI = 20; param.p1 = 0; param.p3 = 10; param.num_load_steps = linspace(0,1,LS); F1crit=403.814571; F1=0.50*F1crit; % problem kinematic and static boundary conditions bc.u_idx = [1 ]; bc.u_val = [0 ]; bc.w_idx = [2 3*N-1]; bc.w_val = [0 0]; bc.phi_idx = [3]; bc.phi_val = [0]; bc.F1_idx = [3*N-2]; bc.F1_val = [-F1]; bc.F3_idx = [ ]; bc.F3_val = [ ]; bc.M_idx = [ ]; bc.M_val = [ ]; % define number of points for evaluating exact solution E = N*100; % obtain exact and approximate solution %[x_exact ,U_exact ] = beam_exact(E, param, bc); [x_exact ,U_exact ] = beam_exact_loadsteps(E, param, bc); [x_approx,U_approx] = beam_approx_loadsteps(N, param, bc); % prepare plotting close all hidden; figure('DefaultAxesFontSize',18); clf; hold on; set(gca,'Ydir','reverse'); % get sub-element values (for plotting and error computation) u_exact = U_exact(1,:); w_exact = U_exact(2,:); u_approx_m = interp1 (x_approx,U_approx(1,:,end), x_exact); w_approx_m = interp1hermite(x_approx,U_approx(2,:,end),U_approx(3,:,end),x_exact); % plot exact and approximate solution graph plot(x_exact ,u_exact , 'b--'); plot(x_exact ,u_approx_m, 'b-' ); plot(x_exact ,w_exact , 'r--'); plot(x_exact ,w_approx_m, 'r-' ); plot(x_approx,U_approx(1,:,end),'bo'); plot(x_approx,U_approx(2,:,end),'ro'); axis([0 param.L -0.005 +0.015]); xlabel('coordinate x [m]'); ylabel('displacement [m]'); title('beam 2nd order theory'); legend('u_{exact}^{2nd order}', 'u_{approx}^{2nd order}', 'w_{exact}^{2nd order}', 'w_{approx}^{2nd order}', 'Location','SouthWest'); grid on; print('beam_2ndorder.png','-dpng'); % plot lambda-deformation curve figure('DefaultAxesFontSize',18); clf; hold on; phi = squeeze(U_approx(3,end,:)); plot(phi, param.num_load_steps, 'b-*'); xlabel('deformation value (end point rotation) [rad]'); ylabel('load factor [-]'); set(gca,'Xdir','reverse'); title('beam 2nd order theory: load-deformation curve'); grid on; print('beam_2ndorder_loadcurve.png','-dpng'); % compute L_2 error w L2_error_w = sqrt( trapz( (w_approx_m - w_exact).^2 ) / trapz(w_exact.^2) ); end
Dernière modification par Antoane ; 19/12/2016 à 22h45. Motif: Ajout balises code