VBA, problema de hidraulica con integrales

Hola,


Hola me han puesto un problema de hidaulica relacionado con integrales y ceros de funciones. Es bastante liante y querria saber si alguien podria ayudarme, tengo en el 1º apartado una duda sobre como estimar el error relativo , y el apartado dos no acavo de programarlo bien:



Calcular la posicion de la seccion cuyo calado es y2 = 2,8 m para un canal trapezoidal
revestido de hormigon (n = 0;014) tal que b = 1;5 m y m = 1, con pendiente I0 = 0,00005 y que transporta un caudal Q = 25,0 m3/s, suponiendo que el calado en la seccion
de referencia es y1 = 2,15 m. Para ello, implementar los metodos del trapecio compuesto
y de Simpson computesto, ya sea como function o como sub. Implementar programas
que consideren discretizaciones del dominio de integracion en puntos equiespaciados, con
n = 2k intervalos, y que determine el valor de k necesario para que el error relativo en
la integral sea menor a una tolerancia jada por el usuario. Dado que el error relativo
no puede calcularse exactamente, emplear las tecnicas de estima del error que vimos en
clase. Para el problema propuesto, >cual es el valor de k si se exige una tolerancia de 10^-6
para los dos metodos empleados?

2. En realidad, en la practica interesa determinar el calado y2 correspondiente a una distancia
L dada, lo que equivale a resolver la ecuacion no lineal
G(y2) = L - integral(de y1 a y2)F(y) dy= 0 (1)
Determinar el calado y2 asociado a una distancia L = 500 m. Para ello implementar los
metodos de biseccion y de Newton. Considerar como tolerancias en los criterios de conver-
gencia 10^-5. Notese que en este caso, es razonable que la tolerancia en la evaluacion de las
integrales sea netamente menor, por ejemplo 10^-7. Mostrar las graficas de convergencia
de ambos metodos. >Que ocurre si al emplear el metodo de Newton la tolerancia en el
calculo de la integral es muy grosera? Comentar los resultados.






Os dejo lo que he programado,para el apartado 1:





Option Explicit

Function func(y As Double) As Double

Dim Q As Double, B_y As Double, G As Double, A_y As Double, I0 As Double, n As Double, Rh As Double, P_y As Double, b As Double, m As Double



Q = Hoja1.Cells(2, 1)

G = Hoja1.Cells(2, 2)

I0 = Hoja1.Cells(2, 3)

n = Hoja1.Cells(2, 4)

m = Hoja1.Cells(2, 5)

b = Hoja1.Cells(2, 6)



B_y = b + 2 * m * y

A_y = (b + m * y) * y

P_y = b + 2 * Sqr(1 + (m ^ 2) * y)

Rh = A_y / P_y


func = ((((Q ^ 2) * B_y) / (G * A_y ^ 3)) - 1) * (I0 - ((n * Q) / (A_y * Rh ^ (2 / 3))) ^ 2) ^ (-1)



End Function


Sub integracion()

Dim metodo As String, y() As Double, n_n As Integer, y1 As Double, y2 As Double, i As Integer, h As Double, integral As Double


metodo = Hoja1.Cells(2, 7)

n_n = Hoja1.Cells(2, 8)

y1 = Hoja1.Cells(2, 9)

y2 = Hoja1.Cells(2, 10)



If metodo = "simpson" And n_n Mod 2 <> 0 Then
MsgBox "Para emplear el método de Simpson el número de intervalos debe ser par"
Cells(5, 2) = "error"
Exit Sub
End If


ReDim y(n_n)
h = (y2 - y1) / n_n
For i = 0 To n_n
x(i) = x1 + i * h
Next i


If metodo = "trapecio" Then
integral = trapecio(x, n)
ElseIf metodo = "simpson" Then
integral = simpson(x, n)
End If
Cells(5, 2) = integral

End Sub



Function trapecio(y() As Double, n_n As Integer) As Double
Dim i As Integer, f1 As Double, f2 As Double
Dim val As Double, h As Double

val = 0


For i = 0 To n_n - 1

f1 = func(y(i))
f2 = func(y(i + 1))
h = (y(i + 1) - y(i))

val = val + (h * (f1 + f2)) / 2

Next i



trapecio = val

End Function



Function simpson(y() As Double, n_n As Integer) As Double

Dim i As Integer, f1 As Double, f2 As Double, f3 As Double
Dim val As Double, h As Double, m_m As Integer


val = 0
m_m = n_n / 2


For i = 1 To m_m

f1 = func(y(2 * i - 2))
f2 = func(y(2 * i - 1))
f3 = func(y(2 * i))

h = y(i + 1) - y(i)

val = val + (h / 3) * (f1 + 4 * f2 + f3)

Next i


simpson = val

End Function





-Para el apartado 2:


Option Explicit


Function func(y As Double) As Double

Dim Q As Double, B_y As Double, G As Double, A_y As Double, I0 As Double, n As Double, Rh As Double, P_y As Double, b As Double, m As Double, L As Double



Q = Hoja1.Cells(2, 1)

G = Hoja1.Cells(2, 2)

I0 = Hoja1.Cells(2, 3)

n = Hoja1.Cells(2, 4)

m = Hoja1.Cells(2, 5)

b = Hoja1.Cells(2, 6)

L = Hoja1.Cells(17, 1)



B_y = b + 2 * m * y

A_y = (b + m * y) * y

P_y = b + 2 * Sqr(1 + (m ^ 2) * y)

Rh = A_y / P_y


func = ((((Q ^ 2) * B_y) / (G * A_y ^ 3)) - 1) * (I0 - ((n * Q) / (A_y * Rh ^ (2 / 3))) ^ 2) ^ (-1)

End Function

Function simpson(y() As Double, n_n As Integer) As Double

Dim i As Integer, f1 As Double, f2 As Double, f3 As Double
Dim val As Double, h As Double, m_m As Integer

n_n = Hoja1.Cells(2, 8)

y1 = Hoja1.Cells(2, 9)

y2 = xk

ReDim y(n_n)
h = (y2 - y1) / n_n
For i = 0 To n_n
x(i) = x1 + i * h
Next i

val = 0
m_m = n_n / 2


For i = 1 To m_m

f1 = func(y(2 * i - 2))
f2 = func(y(2 * i - 1))
f3 = func(y(2 * i))

h = y(i + 1) - y(i)

val = val + (h / 3) * (f1 + 4 * f2 + f3)

Next i


simpson = val

End Function



Sub ceros()

Dim xk As Double, a As Double, tolx As Double, tolf As Double, maxiter As Integer, numiter As Integer
Dim metodo As String

Hoja1.Activate

xk = Cells(2, 9)
tolx = Hoja1.Cells(17, 2)
tolf = Cells(17, 3)
metodo = Cells(17, 4)
maxiter = Cells(17, 6)

If metodo = "Newton" Then
Call MetodoNewton(xk, tolx, tolf, numiter, maxiter)
ElseIf metodo = "Biseccion" Then
a = Cells(17, 5)
Call MetodoBiseccion(xk, a, tolx, tolf, numiter)
Else
MsgBox " Método no reconocido"
numiter = 0

End If

Cells(20, 1) = xk
Cells(20, 2) = numiter
End Sub


Sub MetodoNewton(xk As Double, tolx As Double, tolf As Double, numiter As Integer, maxiter As Integer)

Dim k As Integer, xkmas1 As Double, g_xk As Double, dg_xk As Double, errx As Double, errf As Double, L As Double, val As Double


errx = 2 * tolx
errf = 2 * tolf

k = 0

Cells(20 + k, 10) = k
Cells(20 + k, 11) = xk

Do While (errf > tolf Or errx > tolx) And k <= maxiter

g_xk = L - val
dg_xk = -func(xk)
xkmas1 = xk - g_xk / dg_xk

errx = Abs(xkmas1 - xk) / Abs(xkmas1)
errf = Abs(g_xk)

k = k + 1

xk = xkmas1
Cells(20 + k, 10) = k
Cells(20 + k, 11) = xk
Cells(20 + k, 12) = errx
Cells(20 + k, 13) = errf

Loop

numiter = k
End Sub







Sub MetodoBiseccion(xk As Double, a As Double, tolx As Double, tolf As Double, numiter As Integer)

End Sub