Run Code
|
API
|
Code Wall
|
Misc
|
Feedback
|
Login
|
Theme
|
Privacy
|
Patreon
AjusteNaoLinear-alpha-c
# 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(-5,5,by=0.5);x n=length(x) y=c(76.8, 66.3, 56.4, 47.2, 38.9, 31.4, 24.4, 17.84, 11.1, 3.7, -0.2, 4.1, 11.21, 17.6, 24.1, 31.3, 39.0, 47.33, 56.53, 66.13, 76.73) plot(x,y,col="red") # Por alguma experiência em lidar com dados, supomos ainda que # f_u(s)=(20*s^2+alpha)/(log(as^4+b2*s^2+c)), para algum alpha conhecido. # em que u=(a, b, c) é o vetor de parâmetros desonhecidos. Vamos fixar f(1)=f(x_13)=y_13=1.21 e assim teremos c=e^((20+alpha)/11.21)-a-b # e precisraemos determinar apenas 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_u(x_i))^2, sendp u=(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. alpha=0.1 g<-function(u){ # função g(u) a=u[1];b=u[2];c=exp((20+alpha)/11.21)-a-b p=0 for (i in 1:n){ p=p+( y[i]- (20*x[i]^2+alpha)/log(a*x[i]^4+b*x[i]^2+c) )^2} p/2 } # Teste do código para g(u). g(c(1,2)) gradg<-function(u){ # função gradiente de g(u). a=u[1];b=u[2];c=exp((20+alpha)/11.21)-a-b p=c(0,0) for (i in 1:n){ p=p+( y[i]- (20*x[i]^2+alpha)/log(a*x[i]^4+b*x[i]^2+c))*((20*x[i]^2+alpha)/(log(a*x[i]^4+b*x[i]^2+c))^2 )*( 1/(a*x[i]^4+b*x[i]^2+c) )*c(x[i]^4-1,x[i]^2-1)} p } # Teste da função gradg(u) gradg(c(1,2)) #---- 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*gradg(w)} w } # Teste de aproximação de solução do problema de minimização. u0=c(1,0.5);print("Condição inicial");u0; t=10; m=10000 print("Teste de g(u0)");g(u0) # Teste de g(u0). u=ZeroEuler(u0,t,m) ; print("Aproximação"); u # Aproximadamente u(t) print("Teste de minímo aproximado");g(u) # Teste de minímo aproximado a=u[1];b=u[2];c=exp((20+alpha)/11.21)-a-b print("Aproximação para c"); c # Definição da aproximação de f(u). fap<-function(s){(20*s^2+alpha)/(log(a*s^4+b*s^2+c))} curve(fap,-5,5,ylim=c(0,78)) points(x,y,col="red") # Teste de ajuste. plot(x,y-fap(x),'l') # Se quiser pode aumentar t e n ou escolher outra condição inicial # para se convencer do resultado obtido. # Gráfico para escolha de condição inicial xx <- seq(0.1,2, by=0.05) yy <- seq(1,2.5,by=0.05) require(grDevices) z =matrix(0,length(xx),length(yy)) for ( i in 1:length(xx)){ for (j in 1:length(yy)){ z[i,j]=g(c(xx[i],yy[j]))}} z[is.na(z)] <- 4 op <- par(bg = "white") persp(xx, yy, z, theta = 60, phi = 20, expand = 1) # Gráfico de mf(x,y). contour(xx,yy,z,levels=c(1,10,100,1000))
run
|
edit
|
history
|
help
0
AjusteNaoLinear11-02-21
16-09-2020
K
Julia set variant 3
Example for barplot
EXP1STATS
QRdecomposição
matrix
recursive function
Labegen exemplo