An introduction to the finite element method with applications to nonlinear problems,WHITE,Wiley-Interscience

著書：https://www.amazon.co.jp/s?k=An+introduction+to+the+finite+element+method+with+applications+to+nonlinear+problems&__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&ref=nb_sb_noss

```{r}

#An introduction to the finite element method with applications to nonlinear problems

library(dplyr)

cordi=3

x=sample(1:10,cordi,replace=T)

y=sample(1:10,cordi,replace=T)

u=sample(1:10,cordi,replace=T)

A=cbind(rep(1,cordi),x,y)

alpha=c()

for(j in 1:cordi){
  
A_sub=A

A_sub[,j]=u

alpha=c(alpha,det(A_sub)/det(A))
  
}

N=A%*%alpha/det(A)

u*N


```
```{r}

#p.6

a=0;b=4

t=c(0:4)

m=10;g=9.8

h=0.01

A=array(0,dim=c(length(t)-2,length(t)))

for(j in 1:(length(t)-2)){
  
A[j,j:(j+length(t)-3)]=c(-m/(h^2),(2*m/(h^2)+h),-m/(h^2))  
  
}

#initial condition

vec=c(1,rep(0,length(t)-1))

A=rbind(vec,A)

vec=c(rep(0,length(t)-1),1)

A=rbind(A,vec)

values=c(a,rep(m*g,length(t)-2),b)

y=solve(A)%*%values

```


```{r}

#p.86

library(dplyr)

h=0.01

y=sample(1:10,3)

m=c(0,1,2)

guzai=c(0,0.5,1)*h

A=array(0,dim=c(3,3))

for(j in 1:3){

A[j,]=guzai[j]^m

}

alpha=solve(A)%*%y

X=0

H2=0.01

for(j in 1:length(m)){
  
func=function(x){  
  
f=function(s){
  
return(alpha[j]*x^m[j])  
  
}  

return((((f(x+H)-f(x))/H2)^2+(f(x)^2))/2-32*(f(x)))

}

X=X+integrate(fun,0,h)$value

}


```

```{r}

#p.40

u=sample(1:10,3)

x=sample(1:10,3);y=sample(1:10,3)

A=cbind(rep(1,3),x,y)

alpha=solve(A)%*%u

sum(alpha*c(1,x[1],y[1]))

sum(alpha*c(1,x[2],y[2]))

sum(alpha*c(1,x[3],y[3]))

a1=x[2]*y[3]-y[2]*x[3]

b1=y[2]-y[3];c1=x[3]-x[2]

#N1(x1,y1)==1

sum(c(a1,b1,c1)*c(1,x[2],y[2]))/(det(A))


a2=-x[1]*y[3]+y[1]*x[3]

b2=y[3]-y[1];c2=x[1]-x[3]

#N2(x2,y2)==1

sum(c(a2,b2,c2)*c(1,x[2],y[2]))/det(A)


a3=x[1]*y[2]-y[1]*x[2]

b3=-(y[2]-y[1]);c3=(x[2]-x[1])

#N3(x3,y3)==1

sum(c(a3,b3,c3)*c(1,x[3],y[3]))/det(A)

a=c(a1,a2,a3);b=c(b1,b2,b3);c=c(c1,c2,c3)

#p.48

kij=array(0,dim=c(3,3))

for(i in 1:3){
for(j in 1:3){  
  
kij[i,j]=(b[i]*b[j]+c[i]*c[j])/(4*det(A)/2)  
  
}}

R_e=array(det(A)/(12*2),dim=c(3,3))

diag(R_e)=(det(A)/2)/6

f=sample(1:10,3)

R_e_val=R_e%*%f

```

```{r}

#p.83

A=array(0,dim=c(4,4))

h=0.01

u=c()

for(i in 1:4){
  
vec=sample(1:10,3)  
  
A[i,]=c(1,vec*h)

u=c(u,i)
  
}


d=(svd(A)$d);d[d==0]=10^(-10);d=diag(1/d)

inv_A=svd(A)$u%*%(d)%*%t(svd(A)$v)

alpha=inv_A%*%u

V=det(A)/6

m=1;n=1;l=1;k=1

factorial(m)*factorial(n)*factorial(l)*factorial(k)/(factorial(m+n+l+k+3))*6*V

```

```{r}

#p.90 

x=c(0,2,1);y=c(0,0,1)

r=sqrt((x[2]-x[1])^2+(y[2]-y[1])^2)

a=((x[2]-x[3])*(x[2]-x[1])-(y[3]-y[2])*(y[2]-y[1]))/r

b=((x[3]-x[1])*(x[2]-x[1])+(y[3]-y[1])*(y[2]-y[1]))/r

c=((y[3]-y[2])*(x[2]-x[1])+(x[2]-x[3])*(y[2]-y[1]))/r

mat=array(0,dim=c(6,5))

mat[1,]=c(-b,0,b^2,0,0)
mat[2,]=c((a-b)/2,0,((a-b)^2/4),0,0)
mat[3,]=c(a,0,a^2,0,0)
mat[4,]=c(a/2,c/2,(a/2)^2,a*c/4,(c/2)^2)
mat[5,]=c(0,c,0,0,c^2)
mat[6,]=c(-b/2,c/2,(b^2/4),-b*c/4,(c/2)^2)

A=cbind(rep(1,6),mat)

u=sample(1:10,6)

alpha=solve(A)%*%u

m=c(0,1,0,2,1,0);n=c(0,0,1,0,1,2)

h=function(M_v,N_v){
  
N_v=ifelse(N_v<0,0,N_v);M_v=ifelse(M_v<0,0,M_v)  
  
val=(c^(N_v))*(a^(M_v+1)-(-b)^(M_v+1))*factorial(N_v)*factorial(M_v)/prod(c(1:(N_v+M_v+2)))  
  
return(val)  
   
}

g=array(0,dim=c(6,6))

for(i in 1:6){
for(j in 1:6){  
  
g[i,j]=m[i]*m[j]*h(m[i]+m[j]-2,n[i]+n[j])+n[i]*n[j]*h(m[i]+m[j],n[i]+n[j]-2)  
  
}}

X_u_e=t(alpha)%*%g%*%t(t(alpha))/2

```


```{r}

#p.166

N=10

p=1;C=1;dt=1;K=0.1;h=1/N

a=(p*C*h)/(dt*6)-K/h

b=(p*C*h)/(3*dt)+K/h

u=rep(1,N+1)

mat=array(0,dim=c(N-1,N+1))

for(j in 1:(N-1)){
  
mat[j,j:(j+2)]=c(a,2*b,a)  
  
}

mat=rbind(c(1,rep(0,N)),mat);mat=rbind(mat,c(rep(0,N-1),c(a,b)))

coef=p*C*h/(dt*6)

vec=c(100)

for(j in 1:(N-1)){
  
vec=c(vec,coef*(u[j]+4*u[j+1]+u[j+2]))  
  
}

vec=c(vec,coef*(u[N-1]+2*u[N]))


ite=50

for(s in 1:ite){
  
u_pre=u  
  
u=solve(mat)%*%vec 

vec=c(100)

for(j in 1:(N-1)){
  
vec=c(vec,coef*(u[j]+4*u[j+1]+u[j+2]))  
  
}

vec=c(vec,coef*(u[N-1]+2*u[N]))

print(sum(abs(u-u_pre)))
  
}

```

```{r}

#p.180,181

#simple ordinary differential equation

N=10

A=array(0,dim=c(N-2,N))

for(j in 1:(N-2)){
  
A[j,j:(j+2)]=c(-1,2,-1)  
  
}

A=rbind(c(2,-1,rep(0,N-2)),A)

A=rbind(A,c(rep(0,N-2),-1,2))

f=sample(1:10,N,replace = T)

dx=0.1

A=A/(dx^2)

u=solve(A)%*%f

```

```{r}

#Broyden Quasi-Newton Algorithm p.213

f=function(x){
  
f1=x[1]^2+x[2]^2-1

f2=x[1]^2-x[2]^2+0.5

return(c(f1,f2))
  
}

X=rep(1,2)

ite=100

H=diag(f(X))

for(l in 1:ite){
  
if(sum(abs(f(X)))>10^(-9)){  
  
X_pre=X  
  
X=X-H%*%f(X)  
  
s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

print((f(X)))
  
}
  
}  



```
```{r}

#Broyden Quasi-Newton Algorithm p.213

f=function(x){
  
f1=x[1]^2+x[2]^2-1

f2=x[1]^2-x[2]^2+0.5

return(c(f1,f2))
  
}

X=rep(1,2)

ite=100

H=diag(f(X))

for(l in 1:ite){
  
if(sum(abs(f(X)))>10^(-9)){  
  
X_pre=X  
  
X=X-H%*%f(X)  
  
s=X-X_pre

y=f(X)-f(X_pre)

H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

print((f(X)))
  
}
  
}  



```

```{r}

#pre-Broyden Algorithm p.211

f=function(x){
  
f1=x[1]^2+x[2]^2-1

f2=x[1]^2-x[2]^2+0.5

return(c(f1,f2))
  
}

X=c(1,4)

B=diag((f(X+h)-f(X))/h)

ite=100;h=0.01

for(l in 1:ite){
  
if(sum(abs(f(X)))>10^(-9)){  

X_pre=X

X=X-solve(B)%*%f(X)

s=X-X_pre

y=f(X)-f(X_pre)

B=B+((y-B%*%s)/as.numeric(t(s)%*%s))%*%t(s)

print(f(X))
  
}}


```

```{r}

#Continuation-Newton Algorithm


f=function(x){
  
f1=x[1]^2+x[2]^2-1

f2=x[1]^2-x[2]^2+0.5

return(c(f1,f2))
  
}

X=c(1,4)

ite=100

H=1/ite;h=0.01

for(l in 1:ite){
  
df=(f(X+h)-f(X))/h

X=X-diag(1/df)%*%f(X)

print(f(X))
  
}

```