Sovellus: Differentiaaliyhtälöiden numeeristen ratkaisumenetelmien vertailu

def eulerin_menetelma(f, h, n, t0, y0): f(t, y) = f xy = [(t0, y0)] for k in range(n): tk = xy[-1][0] yk = xy[-1][1] xy.append((tk + h, yk + h*f(tk, yk))) return xy def kp_menetelma(f, h, n, t0, y0): f(t, y) = f xy = [(t0, y0)] for k in range(n): tk = xy[-1][0] yk = xy[-1][1] xy.append((tk + h, yk + h*f(tk + h/2, yk + f(tk, yk)*h/2))) return xy def runge_kutta(f, h, n, t0, y0): f(t, y) = f xy = [(t0, y0)] for k in range(n): tk = xy[-1][0] yk = xy[-1][1] k1 = f(tk, yk) k2 = f(tk + h/2, yk + k1*h/2) k3 = f(tk + h/2, yk + k2*h/2) k4 = f(tk + h, yk + k3*h) k = (k1 + 2*k2 + 2*k3 + k4)/6 xy.append((tk + h, yk + k*h)) return xy @interact def difVertailu(f = input_box(default = "y - 2*sin(t)", label="y' ="), g = input_box(default = "sin(t) + cos(t)", label="Tunnettu ratkaisu: y ="), h = input_box(default=0.5, label="Askelväli: h ="), n = input_box(default=5, label="Askeleita: n ="), t0 = input_box(default=0.0, label="t0 ="), y0 = input_box(default=1.0, label="y0 ="), valinta = selector([(0, "Tunnettu ratkaisu"), (1, "Maxima"), (2, "Ei ratkaista")], label="Tarkkojen arvojen laskeminen"), suunnat=checkbox(default = False, label="Suuntaelementit:")): f(y,t) = f g(y,t) = g # Valitaan, näytetäänkö tarkka ratkaisu (tunnettu tai Maximan määrittämä). if valinta == 0: g(t) = g elif valinta == 1: s = var("s") u = function("u", s) g(t) = desolve(diff(u,s)-f(u,s), u, [t0,y0])(s=t) else: g(t)=0.0 g_kuva = Graphics() if valinta != 2: g_kuva = plot(g(t), (t0, t0+n*h), color="blue", thickness=1.5) # Määritetään ratkaisut eri menetelmien avulla. eu_xy = eulerin_menetelma(f, h, n, t0, y0) kp_xy = kp_menetelma(f, h, n, t0, y0) rk_xy = runge_kutta(f, h, n, t0, y0) # Piirretään ratkaisuun liittyvä grafiikka pisteet = point(eu_xy, color="red") pisteet += point(kp_xy, color="green") pisteet += point(rk_xy, color="brown") viivat = line(eu_xy, color="red") viivat += line(kp_xy, color="green") viivat += line(rk_xy, color="brown") kuva = g_kuva + pisteet + viivat # Suuntaelementtien piirtäminen kuvaajaan. if suunnat == True: kuva += plot_slope_field(f, (t, t0, t0 + n*h), (y, kuva.ymin(), kuva.ymax())) show(kuva) # Tulostetaan taulukko ratkaisufunktioiden arvoista. s = "<br><br><table border=1 cellpadding=5>" s += "<tr><td> $n$ </td><td> $t$ </td><td> $\color{Red}{\\text{Euler}}$ </td><td> $\color{Green}{\\text{Eulerin kp}}$ </td><td> $\color{Brown}{\\text{Runge-Kutta}}$ </td><td> $\color{Blue}{\\text{Tarkka arvo}}$ </td></tr> " for i in range(n+1): s += "<tr><td> $%s$ </td><td> $%s$ </td><td> $%s$ </td><td> $%s$ </td><td> $%s$ </td><td> $%s$ </td></tr>" %(i, N(t0+i*h, digits=4), N(eu_xy[i][1], digits=10), N(kp_xy[i][1], digits=10), N(rk_xy[i][1], digits=10), N(g(t0+i*h), digits=10)) s += "</table>" html(s)