Homepage

Calculus

Not much information on calculus will be logged here, only enough to solve methods like Taylor's.

Derivatives

Integrals

Taylor Series

Runge-Kutta Method

We'll take the equation: a₁ f (t + α₁, y + β₁) and expand it to a 2nd order taylor expansion.

a₁ f(t, y) + a₁ α₁ (THINGf/THINGt) + a₁ β₁ (THINGf/THINGy) + a₁ R₁

T⁽²⁾ = f + (h/2) (THINGf/THINGt) + (h/2) (THINGf/THINGy) f, and a₁ = 1; The coefficients' values are α₁ = (h/2) and β₁ = (h/2) f

After equalizing the coefficients, we have that a₁ = 1, therefore:

T⁽²⁾ = f(t + (h/2), y + (h/2) f(t, y)) - R₁

With this, we achieve a much easier to compute equation, with a recursive function as part of it. The error rate of this equation is of the O(h²) order.

The mid-point Runge-Kutta method (método R-K do ponto médio) is, with ω₀ = α (note, this alpha means initial value, so α = y₀ = y(t₀)):

ω_i+1 = ω_i + h f(t_i + (h/2), ω_i + (h/2) f(t_i, ω_i))

T⁽³⁾ ≈ f(t + α₁, y + δ₁ f(t + α₂, y + δ₂ f(y, t)))

For a fourth order Runge-Kutta equation, things get freaky. Here, we implement a more computable process:

ω₀ = α
K₁ = h f(t_i, ω_i)
K₂ = h f(t_i + (h/2), ω_i + (K₁/2))
K₃ = h f(t_i + (h/2), ω_i + (K₂/2))
K₄ = h f(t_i + h, ω_i + K₃)
ω_i+1 = ω_i + (1/6)(K₁ + 2K₂ + 2K₃ + K₄)

This C code solves 4th order Runge-Kutta equation, from t₀ to t₁, needing only that the user changes the function:

#include
#include

double f(double t, double y)
{
return (log(sin(pow(t, y))) + (pow(t, 2) * pow(y, 2))) / sqrt(t);
}

int main()
{
double w, h, t;

printf("\n**********************\n");
printf("\nMETODO DE RUNGE-KUTTA:\n");
printf("\n**********************\n");

printf("Digite o valor inicial do sistema (alfa): ");
scanf("%lf", &w);

printf("Digite o tempo inicial do sistema (t0): ");
scanf("%lf", &t);

printf("Digite o passo do sistema (h): ");
scanf("%lf", &h);

double K1 = h * f(t, w);
printf("\nK1: %lf", K1);

double K2 = h * f(t + (h/2), w + (K1/2));
printf("\nK2: %lf", K2);

double K3 = h * f(t + (h/2), w + (K2/2));
printf("\nK3: %lf", K3);

double K4 = h * f(t + h, w + K3);
printf("\nK4: %lf", K4);

double W1 = ((K1 + (2*K2) + (2*K3) + K4) * (1.0/6.0)) + w;
printf("\nw1: %lf\n\n", W1);

return 0;
}

Consider a method (R-K or Taylor) whose local truncation error (erro de truncamento local), represtented by τ(h), is of the O(hⁿ) order.

ω₀ = α; ω_i+1 = ω_i + h ϕ(t_i, ω_i , h)

Now consier another method with τ'(h) of O(hⁿ⁺¹) order and approximations:

ω'₀ = α; ω'_i+1 = ω'_i + h ϕ'(t_i, ω'_i , h)

Since ω_i = ω'_i = y(t_i), we can restructure the error function:

τ_i+1(h) = { [ y(t_i+1) - y(t_i) ] / h } - ϕ(t_i, ω_i, h)

This equation then becomes: τ_i+1(h) = [ y(t_i+1) - ω_i+1 ] / h, it being analogous to τ'_i+1(h)

After several transformations, we achieve: τ_i+1(h) = τ'_i+1 + (1/h)(ω'_i+1 - ω_i+1), where τ'_i+1 is a very small value. That is because τ_i+1 is of order N (so, for N = 1, an h of 0.1 would still have 0.1 as product), meanwhile τ'_i+1 is of order N+1 (0.1 to the power of 2 would be 0.01, making it 10 times smaller).

At the end, τ_i+1 ≈ (1/h)(ω'_i+1 - ω_i+1)

In practice, to limit τ_i+1 by ε (a given tolerance), we calculate a value q:

q = [ (ε h) / ( |ω'_i+1 - ω_i+1| ) ]^(1/n);

Note that, here, N = order (So Euler's method would have N = 1, an R-K of order 2 would have N = 2, etc. Euler, R-K and Taylor may be used here).

If q is smaller than 1, then a bad value for h was chosen, likely very big. If q is greater than 1, the entry may be accepted, but the next step should be (q * h). Ideally, q = 1, in which case the best amount of steps (h) was chosen.

Now, obersving the Runge-Kutta-Fehlberg method, it is used with the 4th or 5th order of Runge-Kutta, where each of the Ks is multiplied by specific coefficients. There'd normally be 9 Ks to calculate, but these changes permit only needing to calculate up to K6, with no K2 (although it's still needed to find K3):

ω'_i+1 = ω_i + (16/135)K1 + (6656/12825)K3 + (28561/56430)K4 - (9/50)K5 + (2/55)K6

ω_i+1 = ω_i + (25/16)K1 + (1408/2565)K3 + (2197/4104)K4 - (1/5)K5

The equations for each K are:

holy fuck

Class Notes

These are simply the transcripts (in portuguese) of the class notes my professor has taken, for us to revise as needed.

1. Problema de valor inicial para Equações Diferenciais Ordinárias. A quase totalidade dos fenômenos físicos que observamos ao nosso redor são descritos matematicamente através de Equações Diferenciais. A este procedimento dá-se o nome de modelagem matemática, ou simplesmente modelagem. As Equações Diferenciais resultantes podem ser ordinárias ou parciais, lineares ou não lineares e de ordens diferentes. A grande maioria destas equações não possuem solução analítica e só podem ser resolvidas aproximadamente, utilizando-se métodos numéricos e suas implementações computacionais correspondentes, apropriadas para este fim. No que se segue, iremos estudar métodos de soluções computacionais para EDOs (no tempo), lineares ou não lineares, de primeira ordem. A generalização para sistemas de equações de 1° ordem e para equações de ordem superior será abordada posteriormente.

1.1 - Método de Euler

Apesar do método de Euler ser raramente utilizado na prática, a simplicidade de uma dedução pode ser usada para ilustrar a lógica por trás da dedução das técnicas mais sofisticadas. O objetivo do método é obter uma aproximação para a solução do problema de valor inicial:

(dy/dt) = f(y, t), a ≤ t ≤ b e y(a) = α

Primeiro, vamos investigar graficamente o problema:

Vamos deduzir o método utilizando o teorema de taylor:

y(t_i+1) = y(t_i) + (t_i+1 - t_i) y'(t_i) + [(t_i+1 - t_i)² / 2] y'' (ξ_i)

Para algum ξ_i em [t_i, t_i+1] como h = t_i+1 - t_i

=> y(t_i+1) = y(t_i) + h y'(t_i) + (h/2)² y'' (ξ_i), onde y'(t_i) = f(y, t)

=> y(t_i+1) = y(t_i) + h f(y(t_i), t_i) + (h/2)² y'' (ξ_i), onde y'(t_i) = f(y, t)

Então o método de Euler aproxima y(t_i+1) por excluir termos de ordem superior.

ω₀ = α ; ω_i+1 = ω_i + h f(t_i, ω_i) para cada i = 0, 1, 2, ..., N-1

A equação acima é chamada equação de diferença associada ao método de Euler. Exemplo: y' = y - t² + 1 ; aproxime y para 0 ≤ t ≤ 2 ; y(0) = 0.5. Vamos fazer N = 10, ou seja, h = 0.2, t_i = 0.2 i ; t_N = 2.0 e t₀ = 0.0. Começamos o nosso procedimento computacional fazendo ω₀ = y(t₀) = 0.5. Agora, conforme indicado na equação de diferença associada ao método de Euler: