Lista 5-2023 1

Universidade Federal Rural de Pernambuco - UFRPE UACSA Disciplina: C´alculo Num´erico Professor: Marcelo Flamarion Lista 5 1. Use o M´etodo de Euler para encontrar os primeiros valores y0, y1, y2, y3 para as solu¸c˜oes de cada um dos seguintes problemas de valor inicial. (a) y′ = te3t − 2y, 0 ≤ t ≤ 1, y(0) = 0, com h = 0.1 e h = 0.01. (b) y′ = 1 + (t − y)2, 2 ≤ t ≤ 3, y(2) = 1, com h = 0.25 e h = 0.01. (c) y′ = 1 + y/t, 1 ≤ t ≤ 2, y(1) = 2, com h = 0.1 e h = 0.01. (d) y′ = cos 2t + sin 3t, 0 ≤ t ≤ 1, y(0) = 1, com h = 0.1 e h = 0.01. 2. Repita o exerc´ıcio 1. usando o M´etodo do Trap´ezio. 3. Repita o exerc´ıcio 1. usando o M´etodo de Euler Impl´ıcito. 4. Repita o exerc´ıcio 1. usando o M´etodo de Runge-Kutta de quarta ordem. 5. As solu¸c˜oes dos problemas de valor inicial do exerc´ıcio 1. s˜ao dadas abaixo. Compare os erros absolutos e relativos entre solu¸c˜oes num´ericas e as solu¸c˜oes anal´ıticas. (a) y(t) = 1 5te3t − 1 25e3t + 1 25e−2t. (b) y(t) = t + 1 1−t. (c) y(t) = t ln t + 2t. (d) y(t) = 1 2 sin 2t − 1 3 cos 3t + 4 3. 6. Use o problema teste para determinar a regi˜ao de estabilidade absoluta para o M´etodo do Trap´ezio. import math import numpy as np def f(t,y): return t*math.exp(3*t)-2*y def f_real(t): return t*math.exp(3*t)/5 -math.exp(3*t)/25 + math.exp(-2*t)/25 yi = 0 a = 0 b = 1 h = 0.1 n = int((b-a)/h) t = np.linspace(0,1,n) yreal = [0]*n for i in range (0,n,1): yreal[i] = f_real(t[i]) def euler(yi,a,b,n,t): #questão 1 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi #erroabs[0] = round(abs(y[0] - yreal[0]),5) #errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): y[i] = round(y[i-1] + h*f(t[i-1],y[i-1]),4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("1. A solução pelo método de Euler é: ", y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def trapezio(yi,a,b,n,t): #questão 2 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi #erroabs[0] = round(abs(y[0] - yreal[0]),5) #errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): y[i] = round((y[i-1] + 0.5*h*(f(t[i-1],y[i-1]) + t[i]*math.exp(3*t[i]) ))/(1+h), 4) #muda a cada edo erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("2. A solução pelo método do trapézio é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def eulerimplicito(yi,a,b,n,t): #questão 3 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi #erroabs[0] = round(abs(y[0] - yreal[0]),5) #errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): aux1 = f(t[i-1],y[i-1]) aux2 = f(t[i], y[i-1]+h*aux1) y[i] = round(y[i-1] + (h/2)*(aux1 + aux2),4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("3. A solução pelo método de Euler implícito é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def rk4(yi,a,b,n,t): #questão 4 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi #erroabs[0] = round(abs(y[0] - yreal[0]),5) #errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): k1 = h*f(t[i-1] , y[i-1]) k2 = h*f(t[i-1] + h/2 , y[i-1] + k1/2) k3 = h*f(t[i-1] + h/2 , y[i-1] + k2/2) k4 = h*f(t[i-1] + h , y[i-1] + k3) y[i] = round(y[i-1] + (k1 + k4 + 2*k2 + 2*k3)/6,4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("4. A solução pelo método de RK4 é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") print("t = ", t[0],",",t[1],",",t[2],",",t[3],"\n") euler(yi,a,b,n,t) trapezio(yi,a,b,n,t) eulerimplicito(yi,a,b,n,t) rk4(yi,a,b,n,t) Saída para a) h = 0,1: Para h = 0,01: b) import math import numpy as np from sympy import symbols, solve def f(t,y): return 1 + (t-y)**2 def f_real(t): return t + 1/(1-t) yi = 1 a = 2 b = 3 h = 0.25 n = int((b-a)/h) t = np.linspace(a,b,n+1) yreal = [0]*n for i in range (0,n,1): yreal[i] = f_real(t[i]) def euler(yi,a,b,n,t): #questão 1 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): y[i] = round(y[i-1] + h*f(t[i-1],y[i-1]),4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("1. A solução pelo método de Euler é: ", y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def eulerimplicito(yi,a,b,n,t): #questão 3 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): aux1 = f(t[i-1],y[i-1]) aux2 = f(t[i], y[i-1]+h*aux1) y[i] = round(y[i-1] + (h/2)*(aux1 + aux2),4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("3. A solução pelo método de Euler implícito é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def rk4(yi,a,b,n,t): #questão 4 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): k1 = h*f(t[i-1] , y[i-1]) k2 = h*f(t[i-1] + h/2 , y[i-1] + k1/2) k3 = h*f(t[i-1] + h/2 , y[i-1] + k2/2) k4 = h*f(t[i-1] + h , y[i-1] + k3) y[i] = round(y[i-1] + (k1 + k4 + 2*k2 + 2*k3)/6,4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("4. A solução pelo método de RK4 é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") print("t = ", t[0],",",t[1],",",t[2],",",t[3],"\n") euler(yi,a,b,n,t) #trapezio(yi,a,b,n,t) eulerimplicito(yi,a,b,n,t) rk4(yi,a,b,n,t) OBS o método do trapézio retorna raízes complexas portanto não foi executado Para h = 0,25: H = 0,01: c) import math import numpy as np def f(t,y): return 1 + y/t def f_real(t): return t*math.log(t) + 2*t yi = 2 a = 1 b = 2 h = 0.1 n = int((b-a)/h) t = np.linspace(a,b,n+1) yreal = [0]*n for i in range (0,n,1): yreal[i] = f_real(t[i]) def euler(yi,a,b,n,t): #questão 1 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): y[i] = round(y[i-1] + h*f(t[i-1],y[i-1]),4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("1. A solução pelo método de Euler é: ", y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def trapezio(yi,a,b,n,t): #questão 2 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): aux1 = f(t[i-1],y[i-1]) aux2 = 1 - h/(2*t[i]) y[i] = round((y[i-1] + (h/2)*(1 + aux1))/aux2, 4) #muda a cada edo erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("2. A solução pelo método do trapézio é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def eulerimplicito(yi,a,b,n,t): #questão 3 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): aux1 = f(t[i-1],y[i-1]) aux2 = f(t[i], y[i-1]+h*aux1) y[i] = round(y[i-1] + (h/2)*(aux1 + aux2),4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("3. A solução pelo método de Euler implícito é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def rk4(yi,a,b,n,t): #questão 4 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): k1 = h*f(t[i-1] , y[i-1]) k2 = h*f(t[i-1] + h/2 , y[i-1] + k1/2) k3 = h*f(t[i-1] + h/2 , y[i-1] + k2/2) k4 = h*f(t[i-1] + h , y[i-1] + k3) y[i] = round(y[i-1] + (k1 + k4 + 2*k2 + 2*k3)/6,4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("4. A solução pelo método de RK4 é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") print("t = ", t[0],",",t[1],",",t[2],",",t[3],"\n") euler(yi,a,b,n,t) trapezio(yi,a,b,n,t) eulerimplicito(yi,a,b,n,t) rk4(yi,a,b,n,t) h = 0,1: H = 0,01: d) import math import numpy as np def f(t,y): return math.cos(2*t) + math.sin(3*t) def f_real(t): return 0.5*math.sin(2*t) - (1/3)*math.cos(3*t) + 4/3 yi = 1 a = 0 b = 1 h = 0.1 n = int((b-a)/h) t = np.linspace(a,b,n+1) yreal = [0]*n for i in range (0,n,1): yreal[i] = f_real(t[i]) def euler(yi,a,b,n,t): #questão 1 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): y[i] = round(y[i-1] + h*f(t[i-1],y[i-1]),4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("1. A solução pelo método de Euler é: ", y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def trapezio(yi,a,b,n,t): #questão 2 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): aux1 = f(t[i-1],y[i-1]) aux2 = math.cos(2*t[i]) + math.sin(3*t[i]) y[i] = round((y[i-1] + (h/2)*(aux1+aux2)), 4) #muda a cada edo erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("2. A solução pelo método do trapézio é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def eulerimplicito(yi,a,b,n,t): #questão 3 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): aux1 = f(t[i-1],y[i-1]) aux2 = f(t[i], y[i-1]+h*aux1) y[i] = round(y[i-1] + (h/2)*(aux1 + aux2),4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("3. A solução pelo método de Euler implícito é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") def rk4(yi,a,b,n,t): #questão 4 y = [0]*4 erroabs = [0]*4 errorel = [0]*4 h = (b-a)/n y[0] = yi erroabs[0] = round(abs(y[0] - yreal[0]),5) errorel[0] = round(abs((y[0]-yreal[0])/yreal[0]),5) for i in range (1,4,1): k1 = h*f(t[i-1] , y[i-1]) k2 = h*f(t[i-1] + h/2 , y[i-1] + k1/2) k3 = h*f(t[i-1] + h/2 , y[i-1] + k2/2) k4 = h*f(t[i-1] + h , y[i-1] + k3) y[i] = round(y[i-1] + (k1 + k4 + 2*k2 + 2*k3)/6,4) erroabs[i] = round(abs(y[i] - yreal[i]),5) errorel[i] = round(abs((y[i] - yreal[i])/yreal[i]),5) print("4. A solução pelo método de RK4 é: ",y) print("5. Erro relativo: ", errorel) print("Erro absoluto: ", erroabs,"\n") print("t = ", t[0],",",t[1],",",t[2],",",t[3],"\n") euler(yi,a,b,n,t) trapezio(yi,a,b,n,t) eulerimplicito(yi,a,b,n,t) rk4(yi,a,b,n,t) h = 0,1: H = 0,01: t = 0.0 , 0.01 , 0.02 , 0.03 1. A solução pelo método de Euler é: [1, 1.01, 1.0203, 1.0309] 5. Erro relativo: [0.0, 0.00015, 0.00029, 0.00042] Erro absoluto: [0.0, 0.00015, 0.00029, 0.00043] 2. A solução pelo método do trapézio é: [1, 1.0101, 1.0205, 1.0312] 5. Erro relativo: [0.0, 5e-05, 9e-05, 0.00013] Erro absoluto: [0.0, 5e-05, 9e-05, 0.00013] 3. A solução pelo método de Euler implícito é: [1, 1.0101, 1.0205, 1.0312] 5. Erro relativo: [0.0, 5e-05, 9e-05, 0.00013] Erro absoluto: [0.0, 5e-05, 9e-05, 0.00013] 4. A solução pelo método de RK4 é: [1, 1.0101, 1.0205, 1.0312] 5. Erro relativo: [0.0, 5e-05, 9e-05, 0.00013] Erro absoluto: [0.0, 5e-05, 9e-05, 0.00013] 6.