偏微分方程式の数値シュミレーション　東京大学出版会　登坂宣好・大西和榮先生

```r
#差分解析　p.39
#1次元の数値シュミレーション　
#1次元楕円型方程式の数値シュミレーション
#Strum-Liouville型2点境界値問題


p=function(x){
#return(exp(-x))  
return(1)
}

q=function(x){
#return(x^2)  
return(-1)
}

f=function(x){
#return(exp(-x)*x^2) 
return(0)
}

alpha=1
beta=0

a=0
b=1
n=9
h=(b-a)/(n+1)
X=seq(a,b,h)


mat=array(0,dim=c(n,n))
g=c()
m=c()
d=c()

for(j in 1:n){
  
aj=p(X[j+1]+h/2)+p(X[j+1]-h/2)+q(X[j+1])*h^2
bj=-p(X[j+1]+h/2)

if(j>1 & j<n){

cj=-p(X[j]+h/2)  
gj=f(X[j+1])*h^2  
mat[j,(j-1):(j+1)]=c(cj,aj,bj)
g=c(g,gj)
bj_pre=-p(X[j]+h/2)
mj=cj/dj
m=c(m,mj)
dj=aj-mj*bj_pre
d=c(d,dj)

}else{
  
if(j==n){
  
#j==n:
cj=-p(X[j]+h/2)  
gj=p(X[j+1]-h/2)*beta+f(X[j+1])*h^2    
mat[j,(j-1):(j)]=c(cj,aj)
g=c(g,gj)
bj_pre=-p(X[j]+h/2)
mj=cj/dj
m=c(m,mj)
dj=aj-mj*bj_pre
d=c(d,dj)

}else{
  
#j==1:  
gj=p(h/2)*alpha+f(X[1])*h^2
mat[1,1]=aj
mat[1,2]=bj
g=c(g,gj)
dj=aj
d=c(d,dj)
}}    
}

u=solve(mat)%*%g

v=c(g[1])
for(j in 2:n){
v=c(v,g[j]-m[(j-1)]*v[(j-1)])  
}

```

```r

#p.46 有限要素解析
#1次元の数値シュミレーション　
#1次元楕円型方程式の数値シュミレーション
#Strum-Liouville型2点境界値問題

p=function(x){
#return(exp(-x))  
return(1)
}

q=function(x){
#return(x^2)  
return(1)
}

f=function(x){
#return(exp(-x)*x^2) 
return(sin(x))
}

alpha=1
beta=0

u0=1
a=0
b=1
n=9
h=(b-a)/(n+1)
X=seq(a,b,h)

k=function(I,J){
val=0  
if(I==J){val=(p(X[I])+2*p(X[I+1])+p(X[I+2]))/(2*h)+h*(q(X[I])+6*q(X[I+1])+q(X[I+2]))/12}  
if(J==(I-1)){val=-(p(X[I])+p(X[I+1]))/(2*h)+h*(q(X[I])+q(X[I+1]))/12}
if(J==(I+1)){val=-(p(X[I+1])+p(X[I+2]))/(2*h)+h*(q(X[I+1])+q(X[I+2]))/12}
return(val)
}

mat=array(0,dim=c(n,n))
g=c()
m=c()
d=c()

for(j in 1:n){
  
aj=p(X[j+1]+h/2)+p(X[j+1]-h/2)+q(X[j+1])*h^2

if(j>1 & j<n){
#j>1 & j<n:
mat[j,(j-1):(j+1)]=c(k((j-1),j),k(j,j),k((j+1),j))
gj=h*(f(X[j])+4*f(X[j+1])+f(X[j+2]))/6
g=c(g,gj)
bj_pre=-p(X[j]+h/2)
mj=cj/dj
m=c(m,mj)
dj=aj-mj*bj_pre
d=c(d,dj)

}else{
if(j==n){
#j==n:
mat[j,(j-1):(j)]=c(k((j-1),j),k(j,j))
gj=h*(f(X[j])+4*f(X[j+1])+f(X[j+2]))/6
g=c(g,gj)
bj_pre=-p(X[j]+h/2)
mj=cj/dj
m=c(m,mj)
dj=aj-mj*bj_pre
d=c(d,dj)

}else{
#j==1:  
mat[1,1:2]=c(k(1,1),k(1,2))
gj=h*(f(X[j])+4*f(X[j+1])+f(X[j+2]))/6-(-(p(X[1])+p(X[2]))/(2*h)+h*(q(X[1])+q(X[2]))/12)
g=c(g,gj)
dj=aj
d=c(d,dj)

}}    
}

u=solve(mat)%*%g

v=c(g[1])
for(j in 2:n){
v=c(v,g[j]-m[(j-1)]*v[(j-1)])  
}

```

```r
#p.60 差分解析

a=0
b=1

u0=function(x){
return(sin(pi*x))
}

A=function(tval){
return(0)  
}

B=function(tval){
return(0)  
}

h=0.1
X=seq(a,b,h)

lam=1
dt=0.005
t=0

if(lam*dt/(h^2)<=0.5){

u=c(A(t),u0(X[2:((length(X)-1))]),B(t))

ite=30
for(l in 1:ite){
uvec=c(A(t+dt))
for(j in 2:(length(X)-1)){
uval=u[j]+lam*dt*(u[(j-1)]-2*u[j]+u[(j+1)])/(h^2)
uvec=c(uvec,uval)
}
uvec=c(uvec,B(t+dt))
u=uvec
t=t+dt
plot(X,u)
}

}

```

```r
#p.66 3.2.2 有限要素解析

a=0
b=1

u0=function(x){
return(sin(pi*x))
}

A=function(tval){
return(0)  
}

B=function(tval){
return(0)  
}

h=0.1
X=seq(a,b,h)

lam=1
dt=0.005
t=0

mat1=array(0,dim=c(length(X)-2,length(X)-2))
mat2=array(0,dim=c(length(X)-2,length(X)-2))
for(j in 2:ncol(mat)){
mat1[(j-1),j]=h*(1/6-lam*dt/(h^2))/dt  
mat2[(j-1),j]=1
}
mat1=mat1+t(mat1)
mat2=mat2+t(mat2)
diag(mat1)=h*(4/6+2*lam*dt/(h^2))/dt
diag(mat2)=4

if(lam*dt/(h^2)<=0.5){

u=c(u0(X[2:((length(X)-1))]))

ite=10
for(l in 1:ite){
vector=rep(0,length(u))
vector[1]=A(t)*h*(1/6-lam*dt/(h^2))/dt-(h/(6*dt))*A(t)
vector[length(vector)]=B(t)*h*(1/6-lam*dt/(h^2))/dt-(h/(6*dt))*B(t)
u=solve(mat1)%*%((h/(6*dt))*mat2%*%u-vector)
t=t+dt
plot(X[2:(length(X)-1)],u)
}
}

```

```r
#p.94
#2次元数値シュミレーション
#ポワソン方程式の解法

a=0
b=3
c=4
d=8

Nx=100
Ny=100

hx=(b-a)/Nx
hy=(d-c)/Ny

f=function(X){
return(X[1]^2+X[2]^2)  
}

x=seq(a,b,hx)
y=seq(c,d,hy)
w=0.01

U=function(uvec){
uvec_childer=array(0,dim=c(Nx-1,Ny-1)) 
u11=0;u12=0;u21=0;u22=0
for(i in 1:(Nx-1)){
for(j in 1:(Ny-1)){
if((i+1)==Nx){u11=0}
if((j+1)==Ny){u21=0}
if((j-1)>0){u22=0}
if((i-1)>0){u12=0}
uvec_childer[i,j]=((u11+u12)*hy^2+(u21+u22)*hx^2-f(c(x[i],y[j]))*(hx*hy)^2)/(2*(hx^2+hy^2))
}}
uvec=uvec+w*(uvec_childer-uvec)
return(uvec)
}

u=array(0.1,dim=c(Nx-1,Ny-1))

ite=1000
for(l in 1:ite){
upre=u
u=U(u)  
print(sum(abs(upre-u)))  
}


xvec=c()
yvec=c()
Uvec=c()

x2=x[2:(Nx)]
y2=y[2:(Ny)]

for(i in 1:(Nx-1)){
for(j in 1:(Ny-1)){
xvec=c(xvec,x2[i])
yvec=c(yvec,y2[j])
Uvec=c(Uvec,u[i,j])  
}
}

library(scatterplot3d)
scatterplot3d(xvec, yvec, Uvec)

```