//TP machine n°3 - Résolution numérique d'equa. diff

//Question 1
function z=G(t,y)
    z=3 .*t .*y+1
endfunction

//attention à ne pas appeler "y", le résultat de G. Ici "y" est la variable de G

//Question 2
t=linspace(0,1,100)
a=0 //c'est notre t_0 du cours 
u=0.5 // c'est le y(t_0)
y=ode("rk",u,a,t,G)

fenetre = figure("Figure_name", "Equations", "position", [100 50 1000 600])
fenetre.background = color("white")
set("current_figure", fenetre)
subplot(2,2,1)
//ces quatre dernières lignes sont beaucoup de blabla pour rien
plot2d(t,y,style=[color("black")])

//Question 3

function z=G2(t,y)
    z=3 .*t .*y-1
endfunction

z = ode("rk", u, a, t, G2)
//si je veux la seconde courbe dans une autre fenêtre
subplot(2,2,2)
plot2d(-t, z, style=[color("black")])


//ici, Scilab va superposer dans la fenêtre n°3 les deux courbes. On a donc une courbe sur [-1,1]. J'ai demander de tracer en vert la partie sur [-1,0] et en bleu l'autre.
subplot(2,2,3)
plot2d(t,y,style=[color("blue")])
subplot(2,2,3)
plot2d(-t, z, style=[color("green")])
//pour que ce soit plus joli, on va mettre l'axe de ordonnées au milieu
coor = get("current_axes")
coor.y_location = "origin"

//Question 4
//dans cette question, on va faire varier la condition initiale u
for u=linspace(-1,2,10)//j'ai mis un autre u que dans l'énoncé, c'est plus bô)
    y=ode("rk",u,a,t,G)
    z = ode("rk", u, a, t, G2)
    subplot(2,2,4)
    plot2d(t,y,style=[color("black")])
    plot2d(-t, z, style=[color("black")])
end
coor = get("current_axes")
coor.y_location = "origin"


//Ce qu'il faut retenir de ce TP : définir une fonction (avec la commande "function...endfunction"), utiliser la fonction "ode" , tracer des courbes dans une même fenêtre ou dans des fenêtres différentes (avec "subplot")

//Partie 2 Méthode d'Euler

function [t, y] = euler(g, u, a, b, n)
t = linspace(a, b, n) // calcul du vecteur d’abscisses t(i)
h = (b-a)/(n-1); // c’est la longueur t(i+1)-t(i)
y = zeros(n); // on définit un vecteur contenant n valeurs nulles
y(1) = u; // la première valeur vaut u (condition initiale)
for i = 1:n-1 // cette boucle calcule chaque valeur à partir de la précédente
y(i+1) = y(i) + h * g(t(i), y(i));
end;
endfunction

//Question 5
//on va travailler dans une autre fenêtre graphic pour plus de visibilité
fenetre1 = figure("Figure_name", "Euler", "position", [100 50 1000 600])
fenetre1.background = color("white")
set("current_figure", fenetre1)

[te, ye] = euler(G, 2, 0, 1, 5);
//on a calculé une approximation de la solution sur [0,1], avec la condition initile y(0)=2 et une discretisation de 5 pas. 5 c'est vraiment pas beaucoup, d'où l'aspect "ligne brisée" quand on trace.
//plot2d(te, ye)

//Question 6
[t1, y1] = euler(G, 2, 0, 1, 20);
yref=ode("rk",2,0,linspace(0,1,20),G)
subplot(2,1,1)
plot2d(t1',[y1,yref'], style=[color("green"),color("red")])

//Question 8
N=10:10:500
for i=1:length(N)
    [ta, ya] = euler(G, 2, 0, 1, N(i))
    yex=ode("rk",2,0,linspace(0,1,500),G)
    Err(i)=i*abs(ya(length(ta))-yex(500))
end
subplot(2,1,2)
plot2d(N,Err)
