Run Code
|
API
|
Code Wall
|
Misc
|
Feedback
|
Login
|
Theme
|
Privacy
|
Patreon
Gab1B-newton
# Ajuste de curvas por mínimos quadrados. # Dados os vetores x e y, supomos que y_i=f(x_i)+erro_i, como erro_i um erro "pequeno". x=seq(1930,2010,by=10) x=(x-1930)/10 x # Mudança de escala n=length(x) y=c(0.63,0.76,0.92,1.1,1.2,1.3,1.5,1.7,1.9);y n=length(x) plot(x,y,col="red") # Suponha que y=f(x) (aproximado) e que f(s)= ace^(as)/(1+bce^(as)) e f'(s)=(a-bf(s))f(s). Descobrir aproximações para a, b e c. # # Note que f(0) nos dá c, caso conhecermos a e b. # # Para determinar os valores dos parâmetros resolvemos o problema de minimização # # min g(u) # # em que g(u)= sum (y_i-f(x_i))^2, sendo u=(a,b,c(a,b)). # # Para fazer isso usamos o método de Euler para resolver a equaçõ diferencial # # u'(t)=-grad g(u(t)) # u(0)=u_0 # Seguem os comando para as definições de funções necessárias. g<-function(u){ # função g(u) a=u[1];b=u[2]; c=u[3] p=0 for (i in 1:n){ p=p+( y[i]-a/( b+c*exp(-a*x[i]) ) )^2 } p/2 } # Teste do código para g(u). print(" Teste g"); g(c(1,2,1)) gradng<-function(u){ # função gradiente de g(u). a=u[1];b=u[2]; c=u[3] p=0*u for (i in 1:n){ p=p+( y[i]-a/( b+c*exp(-a*x[i]) ) )*c( -(b+c*exp(-a*x[i]))-a*c*x[i]*exp(-a*x[i]), a, a*exp(-a*x[i]))/(b+c*exp(-a*x[i]) )^2} p } # Teste da função gradg(u) print(" Teste gardiente de g"); gradng(c(1,2,1)) #---- Início da rotin para o método de Euler. ZeroEuler<-function(u0,t,m){ w=u0 h=t/m for (i in 1:m){w=w-h*gradng(w)} w } u0=c(1,1,1) print("Teste condição inicial"); g(u0) u=ZeroEuler(u0,20,20000) #--------------------------------------------------------------------------------------- print("Soulução"); u a=u[1]; b=u[2]; c=u[3] print(" Resíduo da aproximação"); g(u) print("Gradiente no ponto aproximado"); gradng(u) print("Parâmetros a, b e c d encontrados"); a;b;c fap<-function(s){a/(b+c*exp(-a*(s-1930)/10))} # Aproximação da função. print("Gráfico para ajuste") curve(fap,1930,2010) # Teste de ajuste points(seq(1930,2010,by=10),y,col="red") print("População estimada em 2021");fap(2021) # plot(tt,mf(Y[1,],Y[2,]),col = "blue",'l') # curva f(u(t),v(t)) para caminho de maximização obtido pelo método de Euler. #-----------------------Segunda forma: método de Newton para maximizar e resolver a equação grad f(u)=0. hh=0.000001 Hessng<-function(u){ # aproximação de f''(x,y) ou da hessiana de f(x,y). A=matrix(0,3,3) v=u; v[1]=u[1]+hh; w=u;w[1]=u[1]-hh A[1,]=(gradng(v)-gradng(w))/(2*hh) v=u; v[2]=u[2]+hh; w=u;w[2]=u[2]-hh A[2,]=(gradng(v)-gradng(w))/(2*hh) v=u; v[3]=u[3]+hh; w=u;w[3]=u[3]-hh A[3,]=(gradng(v)-gradng(w))/(2*hh) return(A) } Hessng(u) det(Hessng(u)) u=u # condição inicial m=20# número de iterações h=1 # Tamanho do passo for ( i in 1:m){ w=solve(Hessng(u),-h*gradng(u)) u=u+w } print("Aproximação para ponto de mínimo por Newton"); u print("Aproximação para o valor mínimo"); g(u) a=u[1]; b=u[2]; c=u[3] print("Gradiente no ponto aproximado"); gradng(u) print("Hessiana e valores singulares no ponto aproximado");Hessng(u);eigen(Hessng(u))$values; eigen(solve(Hessng(u)))$values print("Parâmetros a, b e c d encontrados"); a;b;c fapb<-function(s){a/(b+c*exp(-a*(s-1930)/10))} # Aproximação da função. print("Gráfico para ajuste") curve(fapb,1930,2010) # Teste de ajuste points(seq(1930,2010,by=10),y,col="red") print("População estimada em 2021");fapb(2021)
run
|
edit
|
history
|
help
0
19-08-2020-GaussSeidelSistema-menos-calc
Gab1A(resumido)
24-08-2020-GradConj
Ch. 5: Cookie Chip Sampling
Call option formula for Laplace distributed outcomes
Tournois Mount & Blade
Critical correlation value calculation
Phyllotaxis in R
3.4.0 / 3.5.2 / 3.5.1 / 3.6.1 / 3.6.5.#
ax² + bx + c = 0 (Solução e gráfico com as raízes)