﻿
#有限要素法で u''+ q(x) u = f(x)  (a < x < b) を解く
#境界条件 u(a) = 1 u(1) = b 

#境界値
global a #グローバル変数として指定
global b
global n
global h
a = 0;
b = 1;
alpha=1;
beta=0;

#a<x<b をn等分
n=10;
h=(b-a)/n;

q = @(x) -1; #q(x)=-1 関数のもう一つの定義法（関数内でもプログラムの変数を引き継ぐ）
f = @(x) -5 * sin(pi*x);


#関数く空間の基底 phi_i(x)
function y=phi(x,i)
  global a #関数内でもグローバル変数として指定
  global n
  global h 
  if (x>= a+(i-1)*h && x <  a+i*h && i>0)
    y=(x- a-(i-1)*h)/h;
  elseif (x>= a+i*h && x <  a+(i+1)*h && i<n)
	y=(a+(i+1)*h-x)/h;
  else
    y=0;
  endif
endfunction

# phi_i(x)の微分
function y=phip(x,i)
  global a
  global n
  global h 
  if (x>= a+(i-1)*h && x <  a+i*h)
    y=1/h;
  elseif (x>= a+i*h && x <  a+(i+1)*h)
	y=-1/h;
  else
    y=0;
  endif
endfunction

#被積分関数
aa=@(x,i,j) -phip(x,i)*phip(x,j)+q(x)*phi(x,i)*phi(x,j);
ff =@(x,i) f(x)*phi(x,i);

#シンプソン積分quad(関数名,積分区間の左端,右端）で計算
A=zeros(n+1,n+1);
for i=0:n
  aaa=@(x) aa(x,i,i);
  A(i+1,i+1)=quad(aaa,0,1);
end
for i=0:n-1
  aaa=@(x) aa(x,i,i+1);
  A(i+1,i+2)=quad(aaa,0,1);
  A(i+2,i+1)=A(i+1,i+2);
end

F=zeros(n+1,1);
for i=1:n-1
 fff=@(x) ff(x,i);
 F(i+1)=quad(fff,0,1);
end

#1行目とn+1行目を整形
F(1)=alpha;
F(n+1)=beta;
A(1,1)=1;
A(1,2)=0;
A(n+1,n)=0;
A(n+1,n+1)=1;

#連立方程式を解く
u=A\F;

plot(linspace(a,b,n+1),u)

#グラフの保存
print -dgif C:\Octave\file.gif
              	






