Gráfica de la solución de un Sistema de Ecuaciones Diferenciales Ordinarias
Publicado por Miguel Angel (1 intervención) el 02/11/2015 06:51:31
Buenas Noches, estoy analizando resolver un sistema de EDO´s y básicamente ya está pero no me compila la gráfica como lo pueden notar estoy tratando de graficar la recta que pasa por los puntos ( j , S ) y ( j+1 , S1 ) esto es en cada iteración pero no me compila. Agradecería su ayuda. Gracias de Antemano.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
Public Class Form1
Dim S, I, R, beta, gamma, S1, I1, R1, S2, I2, R2, S3, I3, R3 As Double
Dim g As Graphics
Function f1(ByVal t As Double, ByVal S As Double, ByVal I As Double, ByVal R As Double) As Double
f1 = -beta * S * I
Exit Function
End Function
Function f2(ByVal t As Double, ByVal S As Double, ByVal I As Double, ByVal R As Double) As Double
f2 = beta * S * I - gamma * I
Exit Function
End Function
Function f3(ByVal t As Double, ByVal S As Double, ByVal I As Double, ByVal R As Double) As Double
f3 = gamma * I
Exit Function
End Function
Private Sub Button1_Click(sender As Object, e As EventArgs) Handles Button1.Click
Dim h, t, k11, k12, k13, k21, k22, k23, k31, k32, k33, k41, k42, k43, k51, k52, k53, k61, k62, k63 As Double
Dim j As Integer
g = PictureBox1.CreateGraphics
g.Clear(Color.White)
g.DrawLine(Pens.Black, 20, 0, 20, 600)
g.DrawLine(Pens.Black, 0, 380, 600, 380)
beta = Val(TextBox1.Text)
gamma = Val(TextBox2.Text)
h = 0.01
S = 999 I = 1 R = 0 For j = 1 To 400 Step 1t = j * h
k11 = h * f1(t, S, I, R)
k12 = h * f2(t, S, I, R)
k13 = h * f3(t, S, I, R)
k21 = h * f1(t + h / 4, S + k11 / 4, I + k12 / 4, R + k13 / 4)
k22 = h * f2(t + h / 4, S + k11 / 4, I + k12 / 4, R + k13 / 4)
k23 = h * f3(t + h / 4, S + k11 / 4, I + k12 / 4, R + k13 / 4)
k31 = h * f1(t + 3 * h / 8, S + ((3 * k11) / 32) + ((9 * k21) / 32), I + ((3 * k12) / 32) + ((9 * k22) / 32), R + ((3 * k13) / 32) + ((9 * k23) / 32))
k32 = h * f2(t + 3 * h / 8, S + ((3 * k11) / 32) + ((9 * k21) / 32), I + ((3 * k12) / 32) + ((9 * k22) / 32), R + ((3 * k13) / 32) + ((9 * k23) / 32))
k33 = h * f3(t + 3 * h / 8, S + ((3 * k11) / 32) + ((9 * k21) / 32), I + ((3 * k12) / 32) + ((9 * k22) / 32), R + ((3 * k13) / 32) + ((9 * k23) / 32))
k41 = h * f1(t + (12 * h / 13), S + ((1932 * k11) / 2197) - ((7200 * k21) / 2197) + ((7296 * k31) / 2197), I + ((1932 * k12) / 2197) - ((7200 * k22) / 2197) + ((7296 * k32) / 2197), R + ((1932 * k13) / 2197) - ((7200 * k23) / 2197) + ((7296 * k33) / 2197))
k42 = h * f2(t + (12 * h / 13), S + ((1932 * k11) / 2197) - ((7200 * k21) / 2197) + ((7296 * k31) / 2197), I + ((1932 * k12) / 2197) - ((7200 * k22) / 2197) + ((7296 * k32) / 2197), R + ((1932 * k13) / 2197) - ((7200 * k23) / 2197) + ((7296 * k33) / 2197))
k43 = h * f3(t + (12 * h / 13), S + ((1932 * k11) / 2197) - ((7200 * k21) / 2197) + ((7296 * k31) / 2197), I + ((1932 * k12) / 2197) - ((7200 * k22) / 2197) + ((7296 * k32) / 2197), R + ((1932 * k13) / 2197) - ((7200 * k23) / 2197) + ((7296 * k33) / 2197))
k51 = h * f1(t + h, S + ((439 * k11) / 216) - (8 * k21) + ((3680 * k31) / 513) - ((845 * k41) / 4104), I + ((439 * k12) / 216) - (8 * k22) + ((3680 * k32) / 513) - ((845 * k42) / 4104), R + ((439 * k13) / 216) - (8 * k23) + ((3680 * k33) / 513) - ((845 * k43) / 4104))
k52 = h * f2(t + h, S + ((439 * k11) / 216) - (8 * k21) + ((3680 * k31) / 513) - ((845 * k41) / 4104), I + ((439 * k12) / 216) - (8 * k22) + ((3680 * k32) / 513) - ((845 * k42) / 4104), R + ((439 * k13) / 216) - (8 * k23) + ((3680 * k33) / 513) - ((845 * k43) / 4104))
k53 = h * f3(t + h, S + ((439 * k11) / 216) - (8 * k21) + ((3680 * k31) / 513) - ((845 * k41) / 4104), I + ((439 * k12) / 216) - (8 * k22) + ((3680 * k32) / 513) - ((845 * k42) / 4104), R + ((439 * k13) / 216) - (8 * k23) + ((3680 * k33) / 513) - ((845 * k43) / 4104))
k61 = h * f1(t + h / 2, S - ((8 * k11) / 27) + (2 * k21) + ((3544 * k31) / 2565) + ((1859 * k41) / 4104) - ((11 * k51) / 40), I - ((8 * k12) / 27) + (2 * k22) - ((3544 * k32) / 2565) + ((1859 * k42) / 4104) - ((11 * k52) / 40), R - ((8 * k13) / 27) + (2 * k23) - ((3544 * k33) / 2565) + ((1859 * k43) / 4104) - ((11 * k53) / 40))
k62 = h * f2(t + h / 2, S - ((8 * k11) / 27) + (2 * k21) + ((3544 * k31) / 2565) + ((1859 * k41) / 4104) - ((11 * k51) / 40), I - ((8 * k12) / 27) + (2 * k22) - ((3544 * k32) / 2565) + ((1859 * k42) / 4104) - ((11 * k52) / 40), R - ((8 * k13) / 27) + (2 * k23) - ((3544 * k33) / 2565) + ((1859 * k43) / 4104) - ((11 * k53) / 40))
k63 = h * f3(t + h / 2, S - ((8 * k11) / 27) + (2 * k21) + ((3544 * k31) / 2565) + ((1859 * k41) / 4104) - ((11 * k51) / 40), I - ((8 * k12) / 27) + (2 * k22) - ((3544 * k32) / 2565) + ((1859 * k42) / 4104) - ((11 * k52) / 40), R - ((8 * k13) / 27) + (2 * k23) - ((3544 * k33) / 2565) + ((1859 * k43) / 4104) - ((11 * k53) / 40))
S1 = S + ((16 * k11) / 135) + ((6656 * k31) / 12815) + ((28561 * k41) / 56430) - ((9 * k51) / 50) + ((2 * k61) / 55)
I1 = I + ((16 * k12) / 135) + ((6656 * k32) / 12815) + ((28561 * k42) / 56430) - ((9 * k52) / 50) + ((2 * k62) / 55)
R1 = R + ((16 * k13) / 135) + ((6656 * k33) / 12815) + ((28561 * k43) / 56430) - ((9 * k53) / 50) + ((2 * k63) / 55)
g.DrawString("Time(días)", Me.Font, Brushes.Black, New Point(300, 383))
g.DrawLine(Pens.Blue, j, S, (j + 1), S1)
'g.DrawLine(Pens.Red, j, I, (j + 1), I1)
'g.DrawLine(Pens.Green, j, R, (j + 1), R1)
S = S1 I = I1 R = R1Next
End Sub
End Class
Valora esta pregunta
0
