Bonjour,

Scilab a un solveur qui s'appelle ODE, êtes-vous autorisé à l'utiliser?
Sinon avec quelle méthode résolution, Euler,Heum?

Bonsoir

Non malheureusement nous n'avons pas le droit d'utiliser ODE.
Nous utilisons la méthode implicite, explicite et des trapèzes ...

Enfaite ce que je voulais savoir et comprendre c'est surtout les paramètres à inscrire dans la fonction f dans le programme ci dessous (une version plus compacte que ce que j'ai mis au dessus) car la résoudre avec les méthodes c'est assez simple mais mon second membre dans mon équation me pose problème sur Scilab...

Merci

Résolution matricielles

OK, en fait la résolution graphique est mauvaise, je ne vois pas nettement votre code. Pourriez-faire un copié-collé afin que je regarde avec Scilab

Oui voila c'est ma résolution qui n'est pas bonne mais je ne comprends pas comment faire pour "intégrer" a scilab le fait que mon second membre n'est pas nul ....
clc
clf
clear

[ok,h]=getvalue("choix du pas",["h"],list("vec",1),["0.05"]);
t=[0:h:50]';

N=prod(size(t))
A=[0,1;-1,-0.5]
B=[0;1]
I=eye(2,2)
f=ones(N,1)
Xexact=zeros(2,N);  Xexpl=zeros(2,N);   Ximpl=zeros(2,N);   Xtrap=zeros(2,N);
    
    for n=2:N
    Xexpl(:,n)=(I+h*A)*Xexpl(:,n-1)+h*B*f(n);
    Ximpl(:,n)=inv(I-h*A)*(Ximpl(:,n-1)+h*B*f(n));
    Xtrap(:,n)=inv(I-h*A/2)*((I+h*A/2)*Xtrap(:,n-1)+h*B/2*(f(n)+f(n-1)));
    Xexact(:,n) = t(n)*(-0.5+t(n)+(exp(-t(n)/4)/2)*cos((sqrt(15)/4)*t(n))-(7/(2*sqrt(15)))*(exp(-t(n)/4)/2)*sin((sqrt(15)/4)*t(n)))+(t(n)-1)*(1-2*(t(n)-1)-(exp(-(t(n)-1)/4))*cos((sqrt(15)/4)*(t(n)-1))+(7/(sqrt(15)))*(exp(-(t(n)-1)/4)/2)*sin((sqrt(15)/4)*(t(n)-1)))+(t(n)-2)*(-0.5+t(n)-2+(exp(-(t(n)-2)/4)/2)*cos((sqrt(15)/4)*(t(n)-2))-(7/(2*sqrt(15)))*(exp(-(t(n)-2)/4)/2)*sin((sqrt(15)/4)*(t(n)-2)))


    end


plot2d(t,[Xexpl(1,',Ximpl(1,',Xtrap(1,',Xexact(1,'],[1,2,3,21])
legends(['Xexpl','Ximpl','Xtrap','Xexact'],[1,2,3,21])
xlabel('t')
ylabel('x(t)')
xgrid();

Je vous met ci joint la résolution


PDF - 421 Ko


PDF - 415 Ko


PDF - 430 Ko

Et voici


PDF - 454 Ko

Ok, merci pour ces éléments, je regarde et reviens vers vous

Pas de soucis encore merci de votre implication

Je viens de comprendre mon erreur, il fallait juste pour définir mon second membre avec une fonction remplie de if...then avec les conditions que j'ai joint ci dessus....
En tout cas merci de m'avoir répondu aussi vite, je pense que vous pouvez supprimer ce topic bien cordialement

Bonjour,

je suis d'accord, j'ai fais tourner votre programme et je l'ai comparer avec le résultat fourni par la fonction ODE, tout est OK sauf Xexact qui a quand même la même allure  à un facteur d'échelle près.
voici votre programme avec la fonction f que j'ai codée

clc
clf
clear
function ft=f(t)
    if t < 1 then
        ft=t // monté
    elseif t >2 then
        ft=0
        else
        ft=2-t // descente
    end
endfunction

[ok,h]=getvalue("choix du pas",["h"],list("vec",1),["0.05"]);
t=[0:h:50]';

N=prod(size(t))
A=[0,1;-1,-0.5]
B=[0;1]
I=eye(2,2)
//f=ones(N,1)
Xexact=zeros(2,N);  Xexpl=zeros(2,N);   Ximpl=zeros(2,N);   Xtrap=zeros(2,N);
    for n=2:N
    Xexpl(:,n)=(I+h*A)*Xexpl(:,n-1)+h*B*f(n*h);
    Ximpl(:,n)=inv(I-h*A)*(Ximpl(:,n-1)+h*B*f(n*h));
    Xtrap(:,n)=inv(I-h*A/2)*((I+h*A/2)*Xtrap(:,n-1)+h*B/2*(f(n*h)+f((n-1)*h)));
    Xexact(:,n) = t(n)*(-0.5+t(n)+(exp(-t(n)/4)/2)*cos((sqrt(15)/4)*t(n))-(7/(2*sqrt(15)))*(exp(-t(n)/4)/2)*sin((sqrt(15)/4)*t(n)))+(t(n)-1)*(1-2*(t(n)-1)-(exp(-(t(n)-1)/4))*cos((sqrt(15)/4)*(t(n)-1))+(7/(sqrt(15)))*(exp(-(t(n)-1)/4)/2)*sin((sqrt(15)/4)*(t(n)-1)))+(t(n)-2)*(-0.5+t(n)-2+(exp(-(t(n)-2)/4)/2)*cos((sqrt(15)/4)*(t(n)-2))-(7/(2*sqrt(15)))*(exp(-(t(n)-2)/4)/2)*sin((sqrt(15)/4)*(t(n)-2)))
    end
n=2:N
length(t)
length(Xexpl(1,:  ))
length(Ximpl(1,:  ))
length(Xtrap(1,:  ))
length(Xexact(1,:  ))
plot2d(t,Xexpl(1,:  ),1)
plot2d(t,Ximpl(1,:  ),2)
plot2d(t,Xtrap(1,:   ),3)
plot2d(t,Xexact(1,:   ),21)
legends(['Xexpl','Ximpl','Xtrap','Xexact'],[1,2,3,21])
xlabel('t')
ylabel('x(t)')
xgrid();

Merci ! J'ai trouvé la meme structure enfaite le facteur d'échelle vient du fait que dans la fonction théorique j'ai mis t(n) t(n-1) mais je voulais faire une fonction heaviside.
Mon code ci vous voulez y jeter un oeil encore merci

clc
clf
clear

Vec_h=[0.05,0.1,0.5]


for i=1:3
    h=Vec_h(i)
    t=[0:h:50]
  

function [f]=fct1(t)    // fonction qui définit le second membre
    [il,ic]=size(t);taille=il*ic;
    f=[];
    for i=1:taille,
        if t(i)<0 then
            f(i)=0;
        end
            if t(i)>0 && t(i)<1 then
                f(i)=t(i);
            end
                if t(i)>1 && t(i)<2 then
                        f(i)=2-t(i);
                end
                     if t(i)>2 then
                         f(i)=0;
                     end
    end
    return
endfunction


function [e]=he1(t)         // fonction d'Heaviside U(T)
    [il,ic]=size(t);taille=il*ic;
    e=[];
    for i=1:taille,
        if t(i)<0 then
            e(i)=0;
        end
        if t(i)>0 then
            e(i)=1;
        end
        end
    return
            
endfunction

function [g]=he2(t)         // fonction d'Heaviside U(T-1)
    [il,ic]=size(t);taille=il*ic;
    g=[];
    for i=1:taille,
        if t(i)<1 then
            g(i)=0;
        end
        if t(i)>1 then
            g(i)=1;
        end
    end
    return
            
endfunction

function [j]=he3(t)         // fonction d'Heaviside U(T-2)
    [il,ic]=size(t);taille=il*ic;
    j=[];
    for i=1:taille,
        if t(i)<2 then
            j(i)=0;
        end
        if t(i)>2 then
            j(i)=1;
        end
        end
    return
            
endfunction

//définition de parametres pour alléger la notation
f=fct1(t)
a=he1(t)
b=he2(t)
c=he3(t)

N=prod(size(t))
A=[0,1;-1,-0.5]
B=[0;1]
I=eye(2,2)

Xexact=zeros(2,N);  Xexpl=zeros(2,N);   Ximpl=zeros(2,N);   Xtrap=zeros(2,N);
    
    for n=2:N
        
    Xexpl(:,n)=(I+h*A)*Xexpl(:,n-1)+h*B*f(n);       // fonction explicite
    
    Ximpl(:,n)=inv(I-h*A)*(Ximpl(:,n-1)+h*B*f(n));  // fonction implicite
    
    Xtrap(:,n)=inv(I-h*A/2)*((I+h*A/2)*Xtrap(:,n-1)+h*B/2*(f(n)+f(n-1))); //fonction trapèzoidale
    
    Xexact(:,n) = a(n)*(-0.5+t(n)+(exp(-t(n)/4)/2)*cos((sqrt(15)/4)*t(n))-(7/(2*sqrt(15)))*(exp(-t(n)/4))*sin((sqrt(15)/4)*t(n)))+b(n)*(1-2*(t(n)-1)-(exp(-(t(n)-1)/4))*cos((sqrt(15)/4)*(t(n)-1))+(7/(sqrt(15)))*(exp(-(t(n)-1)/4))*sin((sqrt(15)/4)*(t(n)-1)))+c(n)*(-0.5+t(n)-2+(exp(-(t(n)-2)/4)/2)*cos((sqrt(15)/4)*(t(n)-2))-(7/(2*sqrt(15)))*(exp(-(t(n)-2)/4))*sin((sqrt(15)/4)*(t(n)-2)))                                  // fonction exacte

    end

subplot(3,3,3+i)
plot2d(t,[Xexpl(1,',Ximpl(1,',Xtrap(1,',Xexact(1,'],[1,2,3,21])
legends(['Xexpl','Ximpl','Xtrap','Xexact'],[1,2,3,21])
xlabel('t')
ylabel('x(t)')
xgrid();
end

salut,
pour echanger du code il y a sous le message le bouton </>

@alb12 ,ok merci beaucoup
@Mars1933, je vais essayer votre programme dans Scilab

je pense qu'il reste encore un souci pour  Xexpl, j'ai un doute aussi sur les autres à causes des valeurs.
voici 2 images, une de votre programmes avec les 4 courbes, l'autre de la version que je vous ai passée sans Xexact

Enfaite le probleme vient du fait que vous ne déclarez pas de fonction d'Heaviside comme j'ai fais l'erreur. Au lieu de mettre t(n) t(n-1) etc il faudrait mettre U(n) U(n-1) comme le programme que je vous ai renvoyé.


function [e]=he1(t)         // fonction d'Heaviside U(T)
    [il,ic]=size(t);taille=il*ic;
    e=[];
    for i=1:taille,
        if t(i)<0 then
            e(i)=0;
        end
        if t(i)>0 then
            e(i)=1;
        end
        end
    return
            
endfunction

function [g]=he2(t)         // fonction d'Heaviside U(T-1)
    [il,ic]=size(t);taille=il*ic;
    g=[];
    for i=1:taille,
        if t(i)<1 then
            g(i)=0;
        end
        if t(i)>1 then
            g(i)=1;
        end
    end
    return
            
endfunction

function [j]=he3(t)         // fonction d'Heaviside U(T-2)
    [il,ic]=size(t);taille=il*ic;
    j=[];
    for i=1:taille,
        if t(i)<2 then
            j(i)=0;
        end
        if t(i)>2 then
            j(i)=1;
        end
        end
    return
            
endfunction


On voit en zoomant sur chaque tracé que la courbe qui s'approche le plus de celle théorique et celle avec l'approximation trapézoïdale celon le pas h cela varie et la méthode explicite et implicite devienne juste incohérente

Je me trompe peut-être mais j'ai un doute sur ce vous trouvez pour Xexpl

 for n=2:N
    Xexpl(:,n)=(I+h*A)*Xexpl(:,n-1)+h*B*f((n-1)*h);
    Ximpl(:,n)=inv(I-h*A)*(Ximpl(:,n-1)+h*B*f((n-1)*h));
    Xtrap(:,n)=inv(I-h*A/2)*((I+h*A/2)*Xtrap(:,n-1)+h*B/2*(f((n-1)*h)+f((n-2)*h)));

erratum lire : ce n'est pas vrai dans votre code.

Quand vous démarré la simulation on vous demande la valeur de h, dans le cadre de mon exercice on me demande h=0,05 puis 0,1 puis 0,5 et en faisant la simulation avec votre programme qui est totalement juste mais en intégrant que h=0,5 on retrouve bien les memes resultats, avec cette oscillation de Xexpl

Ok,ça marche, je n'avais pas vu qu'on vous demandait plusieurs valeur pour h.