[MatLab] Résolution d'équations différentielles et autres problèmes

[invitede4c85a1]

Bonjour,

Je suis plus que débutant en Matlab, et à vrai dire j'ai toujours été allergique à ce genre de logiciels (Matlab, Maple, etc...)
J'ai récemment du m'initier à LabView et ça ne m'a posé de soucis par contre, même si le but final diffère.
En tous cas, je suis à la recherche d'une aide pour comprendre certains programmes. J'ai essayé de comprendre leur but et comment ils fonctionnaient mais je suis complétement paumé. Je ne peux pas utiliser MatLab chez moi, n'ayant déjà pas la Symbolic ToolBox.

Mon but est de réussir à résoudre ces problèmes.

#### liens supprimés

J'ai mis en pièces jointes les différents programmes à ma disposition.
Pour le problème 1, faut-il utiliser les itérations newtoniennes et pour le problème 2 le programme Euler ?

Merci d'avance.

Itérations de Newton : #### liens supprimés
Gauss : #### liens supprimés
Methode des moindres carrés Symbolic : #### liens supprimés
Euler : #### liens supprimés
Méthode des moindres carrés Numerical : #### liens supprimés
-----

[Jack]

Je ne peux pas utiliser MatLab chez moi, n'ayant déjà pas la Symbolic ToolBox.

tu as essayé avec octave, le clone de matlab?

[invitede4c85a1]

J'ai essayé d'installer Octave mais je n'arrive pas à le lancer.
Je ne peux juste pas utiliser les variables symboliques sous Matlab.

J'ai réussi le a) de mon exercice 1 grâce à la méthode du pivot de Gauss.
Mais comment faire pour le b) vu la présence des puissances ?

J'aimerais bien une petite indication pour l'exercice 2 également, merci.

[JPL]

Les images doivent être postées en pièces jointes et les codes doivent être postés dans le corps des messages en utilisant la balise Code (# dan la barre d'outils du mode avancé). Merci.

[A voir en vidéo sur Futura]

[invitede4c85a1]

EULER

%  Pondilium equation x - angle, x=0 - minimum of potential energy 
%  d^2(x)/ dt^2 + w0^2 * sin(x) = 0 
%       =>  
%  [ dx/dt = y; dy/dt = - sin(x) ] y - angular speed
%       =>
%  dq/dt = f(q); q = [ x; y ]  f(q) = [y; -sin(x)]

f = @(q) [ q(2); -sin(q(1)) ];  % function to integrate, 
                                % it could also be a function of t 

% first experiment, small step - the both OK
%q0 = [0; 1.5];  % initial point: oscilations
                 % min of potential energy, max of kinetical energy
%dt = 0.001; % timestep
%steps = 10000; % number of iterations

% second experiment, small step - the both OK
q0 = [0; 3];  % initial point: pondilium is turning...
              % min of potential energy, max of kinetical energy
dt = 0.001; % timestep
steps = 10000; % number of iterations

% third experiment - Euler divergies after one turn
% q0 = [0; 1.5];  % initial point
% dt = 0.2; % timestep
% steps = 100; % number of iterations

% forth experiment - Both diverges
% q0 = [0; 1.5];  % initial point
% dt = 1; % timestep
% steps = 20; % number of iterations

% simple Euler method
q_sEuler = NaN(size(q0,1), steps);
q_sEuler(:,1)=q0;

for i=1:steps-1
  q_sEuler(:,i+1) =  q_sEuler(:,i) + f(q_sEuler(:,i))*dt; 
end;

plot(q_sEuler(1,:), q_sEuler(2,:), 'b');
hold on;

% modified Euler method
q_mEuler = NaN(size(q0,1), steps);
q_mEuler(:,1)=q0;

for i=1:steps-1
  t_i_plus05 = (i-1) * dt + 1/2 * dt;
  q_i_plus05 = q_mEuler(:,i) + 1/2 * f(q_mEuler(:,i))*dt;
  f_i_plus05 = f(q_i_plus05);
  q_mEuler(:,i+1) =  q_mEuler(:,i) + f_i_plus05*dt; 
end;

plot(q_mEuler(1,:), q_mEuler(2,:), 'r');

LSF Numerical

me=5.485799E-4;
Z=[1;1;1;2;2;3;3;4;5;5;6;6;7;7;8;8;8;13;15;19;19;21;22;27;35;39;45;55;57;65;69;70;73;79;83;92;94];
np=[1;2;3;3;4;6;7;9;10;11;12;13;14;15;16;17;18;27;31;39;41;45;48;59;79;89;103;133;139;159;169;171;181;197;209;238;238];
m=[1.007825;2.01410;3.01605;3.01603;4.00260;6.01512;7.01600;9.01218;10.0129;11.093;12.000;13.00335;14.00307;15.00011;15.99491;16.99913;17.99917;26.9815;30.97376;38.9637;40.9618;44.95591;47.94795;58.93320;78.9183;88.90586;102.90550;132.90543;138.90636;158.92535;168.93423;170.9363;180.94801;196.96656;208.98039;238.05079;238.04956];

Meffexp=(m-me*Z)./np;
% ./ est l'opérateur de division vectorielle

% fitting function f(Z)=a+b*Z+c*Z^2+d/(Z^Kval) , Kval=0.5

N=size(Z,1);
y=zeros(4,1);
for i=1:N
    y(1)=y(1)+2*Meffexp(i);
    y(2)=y(2)+2*Meffexp(i)*Z(i);
    y(3)=y(3)+2*Meffexp(i)*(Z(i))^2;
    y(4)=y(4)+2*Meffexp(i)*(Z(i))^(-0.5);
end;
A=zeros(4,4);
for i=1:N;
    A(1,1)=A(1,1)+2;
    A(1,2)=A(1,2)+2*Z(i);
    A(1,3)=A(1,3)+2*(Z(i))^2;
    A(1,4)=A(1,4)+2*(Z(i))^(-0.5);
    A(2,1)=A(2,1)+2*Z(i);
    A(2,2)=A(2,2)+2*Z(i)^2;
    A(2,3)=A(2,3)+2*(Z(i))^3;
    A(2,4)=A(2,4)+2*(Z(i))^(0.5);
    A(3,1)=A(3,1)+2*(Z(i))^2;
    A(3,2)=A(3,2)+2*(Z(i))^3;
    A(3,3)=A(3,3)+2*(Z(i))^(4);
    A(3,4)=A(3,4)+2*(Z(i))^(1.5);
    A(4,1)=A(4,1)+2*(Z(i))^(-0.5);
    A(4,2)=A(4,2)+2*(Z(i))^(0.5);
    A(4,3)=A(4,3)+2*(Z(i))^(1.5);
    A(4,4)=A(4,4)+2*(Z(i))^(-1);
end

abcd=mygausspivot(A,y);

Mefffit=abcd(1)+abcd(2)*Z+abcd(3)*(Z.*Z)+abcd(4)./sqrt(Z);

plot_handle = plot(Z, Meffexp,' ob', Z, Mefffit,'-r');

set(plot_handle(1),'DisplayName','Experiment data');
set(plot_handle(2),'DisplayName','Fitting data');
legend('show');
xlabel('Charge');
ylabel('Effective Mass');
title('Effectve mass and nuclear stability') 

disp('Mean square deviation:')
disp(sqrt( (Meffexp-Mefffit)' * (Meffexp-Mefffit) )/size(Z,1));

LSF SYMBOLIC

me=5.485799E-4; % mass of electron in atomic mass units

% elements data
Z=[1;1;1;2;2;3;3;4;5;5;6;6;7;7;8;8;8;13;15;19;19;21;22;27;35;39;45;55;57;65;69;70;73;79;83;92;94];
np=[1;2;3;3;4;6;7;9;10;11;12;13;14;15;16;17;18;27;31;39;41;45;48;59;79;89;103;133;139;159;169;171;181;197;209;238;238];
m=[1.007825;2.01410;3.01605;3.01603;4.00260;6.01512;7.01600;9.01218;10.0129;11.093;12.000;13.00335;14.00307;15.00011;15.99491;16.99913;17.99917;26.9815;30.97376;38.9637;40.9618;44.95591;47.94795;58.93320;78.9183;88.90586;102.90550;132.90543;138.90636;158.92535;168.93423;170.9363;180.94801;196.96656;208.98039;238.05079;238.04956];

Meffexp=(m-me*Z)./np;  % effective mass calculation

%Meff = @(Z, n, m) (m-me*Z)/n;      % online effective mass function

Kval = 0.1; % fix Kval parameter

syms a b c d;
Feff = @(Z) a+b*Z+c*Z^2+d/Z^Kval;   % online fiting function
%syms S(a, b, c, d);    % define S as a function of a,b,c,d - MATLAB >= 2011
syms S; % MATLAB 2010
S = 0;
for i=1:size(Z,1)
    S = S + (Feff(Z(i))-Meffexp(i))^2;
end;
S = simplify(vpa(S));  % vpa calcuates the fractions, powers, etc...
%dS = gradient(S);      % dS is a set of 4 functions  MATLAB >= 2011
dS = [diff(S, a); diff(S, b); diff(S, c); diff(S, d)]; % MATLAB 2010
ss = solve(dS);        % find the solution dS==0, the same thing that :
                       % ss = solve(dS(1)==0, dS(2)==0,  dS(3)==0,  dS(4)==0, a, b, c, d )

disp('Retrieved parameters:' );
disp([ss.a ss.b ss.c ss.d]); % print out the parameters

Mefffit = ss.a + ss.b*Z + ss.c*(Z.*Z)+ ss.d ./ sqrt(Z);

plot_handle = plot(Z, Meffexp,' ob', Z, Mefffit,'-r');

set(plot_handle(1),'DisplayName','Experiment data');
set(plot_handle(2),'DisplayName','Fitting data');
legend('show');
xlabel('Charge');
ylabel('Effective Mass');
title('Effectve mass and nuclear stability') 

disp('Mean square deviation:')
disp(double(sqrt(subs(S, symvar(S), [ss.a ss.b ss.c ss.d]) )/size(Z,1)));

GAUSS

function X=mygauss(A,b)
    n=size(A,1);
    X = NaN(n,1);
    
    % Matrix triangualisation
    for i=1:n
        
        p=A(i,i); % pivoting element
        for k=i:n
            A(i,k)=A(i,k)/p;
        end
        b(i)=b(i)/p;        
        for j=i+1:n
            r=A(j,i);
            for k=i:n
                A(j,k)=A(j,k)-A(i,k)*r;
            end
            b(j)=b(j)-b(i)*r;
        end
    end
    
    % X calculation - triangular matrix to unit matrix
    for i=n-1:-1:1
        for k=i+1:n
            b(i)=b(i)-b(k)*A(i,k);
        end
    end
    X=b;

end

NEWTONIAN ITERATIONS

% Resolution of non-liniar system of equations 
% by the newtonean iterations
clear;

% First example +
syms x1 x2;    		% define two symbolic variables x1 and x2.
xi = [x1; x2];         % vector of variables
syms f(x1, x2)  	% define the simbolic fonction. Works on MATLAB >= 2011
syms f; 		% OK for MATLAB <= 2010
f = [x1^3-3*x1*x2^2-1; 3*x1^2*x2-x2^3 ]; % function f: R2 -> R2
%f = [x1+x2-2; x1-x2] % other examples
%f = [x1^2+x2^2-8;x1^2-x2^2]
xi0 = [-0.5; 0.5];      % initial point first solution
%xi0 = [-0.5; -0.5];      % initial point second solution
%xi0 = [0.6; -0.01];      % initial point theard solution

% Second and theard examples
%syms x1 x2 x3;    	 % define a symbolic variables.
%xi = [x1; x2; x3];      % variables
%f=[x1^2+x2^2+x3^2-12; x1+x2-4; x1+x3-4];
%%f=[x1^2+x2^2+x3^2-48;x1+x2-8;4*x1+4*x2-20] % error in J invertion as second and theurd equations are incompartable 
%xi0 = [-0.5; 0.5; 0.5]; % initial point

% % Forth example
% syms xi;
% f = xi^3-2*xi+2;
% xi0 = 0;              % infinite loop for this initial value

J=jacobian(f);          % Jacobi matrix of f

% iJ = inv(J);            % inverse Jacobian matrix
% df = iJ * f;            % correction vector f(x)/f’(x)
df = J\f;               % better implementation that two previouse lines

f0 = subs(f, xi, xi0);  % function value at initial point
nf0 = norm(f0, 2);
tau = 10^-2;            % relative tolerance of solution improvement
eps = 10^-2;            % xi - x(i-1) tolerance
imax = 32;              % maximum number of iterations

disp('i                  xi                 norm(xi-x(i-1))                f(xi)                     norm(f(xi)) ');
disp( ['0            [',  num2str(xi0'), ']              none            [', num2str(f0'), ']       ', num2str(nf0) ]);

for i=1:imax
  xi1 = xi0 - subs(df, xi, xi0);  % x1 = x0 - J^(-1)*f
  f1 = subs(f, xi, xi1);
  nf1 = norm(f1, 2);
  disp( [ num2str(i), '           [',  num2str(xi1'), ']           ', num2str(norm(xi0-xi1)),'           [', ...
      num2str(f1'), ']       ', num2str(nf1) ]);
  if (nf1==0) || ( abs(nf0/nf1-1) < tau ) % exit loop if norm(f) does not become smaller //// signification de ||  :  "ou"
      disp('Exit by tau or by zero norm...');
      break;      
  end;
  if norm(xi0-xi1)<eps
      disp('Exit by eps...');
      break;
  end
  xi0=xi1;
  nf0=nf1;
end;

[invitede4c85a1]

Je pense avoir réussi l'exercice 2 en partie.

J'ai utilisé une solution commune de l'équation différentielle, que je devrais normalement obtenir à partir d'un programme Matlab.

a) L'énergie potentielle est maximale si la hauteur est maximale, soit en x = -1m.
b) L'énergie cinétique est maximale si l'énergie potentielle est minimale, soit en x = 1m.

J'ai mis en pièce jointe le schéma permettant de comprendre et la figure de mouvement en phase (x,dx/dt).

Quel valeur de pas dois-je choisir pour le temps pour avoir une petite erreur sans pour autant trop en demander à Matlab ?