restart;
phi:=proc(x,L)
if x>L then
0
else
if x>=0 then 1-x else
phi(-x,L)
end if;
end if;
end proc;
phiLocal:=proc(i,x,h,L)
phi(x/h-(i-1),L)
end proc;
yapprox:=proc(x,c,h,n,L)
local i;
add(c[i]*phiLocal(i,x,h,L),i=1..n)
end proc;
L:=1; nPoints:=100; nElements:=nPoints-1;
q:=4; f:=4;
nShapeFunctions:=nPoints;
h:=evalf(L/nElements);
#plot([''phiLocal(i,x,h,L)'' $i=1..nPoints],x=0..L);
with(LinearAlgebra):
makeMatrix:=proc(q,h,f,nPoints)
local A,i,j;
A:=Matrix(nPoints);
for i from 1 to nPoints do
for j from 1 to nPoints do
if i=j then
A[i,j]:=2+2/3 * q*h^2;
else
if j=i-1 or j=i+1 then
A[i,j]:= -1+1/6 *q *h^2;
end if;
end if;
end do;
end do;
A
end proc:
A:=makeMatrix(q,h,f,nPoints):
b:=Vector(nPoints):
b[1..-1]:=h^2*f:
c:=LinearSolve(A,b):
with(plots):
listplot(['yapprox(i*h,c,h,nPoints,L)' $i=0..nPoints-1]);
Zio2JEkieEc2IkkiTEdGJUYlRiVGJUAlMjklOSQiIiFAJTFGK0YqLCYiIiJGL0YqISIiLUkkcGhpR0YlNiQsJEYqRjBGKUYlRiVGJQ==
Zio2JkkiaUc2IkkieEdGJUkiaEdGJUkiTEdGJUYlRiVGJS1JJHBoaUdGJTYkLCgqJjklIiIiOSYhIiJGLzkkRjFGL0YvOSdGJUYlRiU=
Zio2J0kieEc2IkkiY0dGJUkiaEdGJUkibkdGJUkiTEdGJTYjSSJpR0YlRiVGJS1JJGFkZEclKnByb3RlY3RlZEc2JComJjklNiM4JCIiIi1JKXBoaUxvY2FsR0YlNiZGNDkkOSY5KEY1L0Y0O0Y1OSdGJUYlRiU=
IiIi
IiQrIg==
IiMqKg==
IiIl
IiIl
IiQrIg==
JCIrNTU1NTUhIzY=
NigtJSdDVVJWRVNHNiM3YHE3JCQiIzUhIiIkIjIjW0o9PDc+TjohIz03JCQiIz8hIiIkIjJDWlFNMyU+SUkhIz03JCQiI0khIiIkIjEkb0tLeT1jWyUhIzw3JCQiI1MhIiIkIjEkM0hzTWY/IWYhIzw3JCQiI10hIiIkIjEvM3UlKlE0IUcoISM8NyQkIiNnISIiJCIxTykqb3BbRz8nKSEjPDckJCIjcSEiIiQiMWMibyxGekojKiohIzw3JCQiIyEpISIiJCIyXFlwcikzJCo9NiEjPDckJCIjISohIiIkIjJbKmVVVyE+PkMiISM8NyQkIiQrIiEiIiQiMnMyQ2pmSzhPIiEjPDckJCIkNSIhIiIkIjJlM0g5Rz9zWiIhIzw3JCQiJD8iISIiJCIyOHJEIypSSCcqZSIhIzw3JCQiJEkiISIiJCIyYnJCJUhlZylwIiEjPDckJCIkUyIhIiIkIjJZZyRvXVM+Lz0hIzw3JCQiJF0iISIiJCIydmNxKGVyVjE+ISM8NyQkIiRnIiEiIiQiMk0pUjklKW9QMD8hIzw3JCQiJHEiISIiJCIxJW8qKSozTzAsQCEjOzckJCIkIT0hIiIkIjJjIj5xJFExTj4jISM8NyQkIiQhPiEiIiQiMkd0QEclSHgjRyMhIzw3JCQiJCsjISIiJCIxRyRRL3MqKSlvQiEjOzckJCIkNSMhIiIkIjIzZWZdJz0qPVgjISM8NyQkIiQ/IyEiIiQiMicpKnkmUkQ4PWAjISM8NyQkIiRJIyEiIiQiMnc/KCpwXSdvM0UhIzw3JCQiJFMjISIiJCIxQm4lKioqSGEjbyMhIzs3JCQiJF0jISIiJCIyYSZwR3hHVGBGISM8NyQkIiRnIyEiIiQiMm06OFgxRDgjRyEjPDckJCIkcSMhIiIkIjJsQ1QqenNJJylHISM8NyQkIiQhRyEiIiQiMUA1LVlnUVtIISM7NyQkIiQhSCEiIiQiMXg3OytuZTJJISM7NyQkIiQrJCEiIiQiMWU3MjBNJFIxJCEjOzckJCIkNSQhIiIkIjFcRWllIlx1NiQhIzs3JCQiJD8kISIiJCIxVmhCLmU6b0ohIzs3JCQiJEkkISIiJCIwVCd6TVMyO0shIzo3JCQiJFMkISIiJCIxNiIqNDZNQWhLISM7NyQkIiRdJCEiIiQiMicpUkwpZkJpLkwhIzw3JCQiJGckISIiJCIybTsuaD0pR1ZMISM8NyQkIiRxJCEiIiQiMT0sXHpxQiFRJCEjOzckJCIkIVEhIiIkIjFLOm0/VFs5TSEjOzckJCIkIVIhIiIkIjFFYV8oR1ZnVyQhIzs3JCQiJCslISIiJCIxT2MkNFlGXFokISM7NyQkIiQ1JSEiIiQiMUMiXCpIJVs2XSQhIzs3JCQiJD8lISIiJCIxeHdqJypvckNOISM7NyQkIiRJJSEiIiQiMmIlZlghW1VjYSQhIzw3JCQiJFMlISIiJCIxQ3c7QVAkUmMkISM7NyQkIiRdJSEiIiQiMWM4TCgzKWZ6TiEjOzckJCIkZyUhIiIkIjJkelokcD5rI2YkISM8NyQkIiRxJSEiIiQiMXQhcD9wcUlnJCEjOzckJCIkIVshIiIkIjJ4cnI+XikpM2gkISM8NyQkIiQhXCEiIiQiMngkUiopPicpNDtPISM8NyQkIiQrJiEiIiQiMCd6TFVKcT1PISM6NyQkIiQ1JiEiIiQiMjEnekxVSnE9TyEjPDckJCIkPyYhIiIkIjJ4JFIqKT4nKTQ7TyEjPDckJCIkSSYhIiIkIjE8PCg+XikpM2gkISM7NyQkIiRTJiEiIiQiMjkycD9wcUlnJCEjPDckJCIkXSYhIiIkIjJCelokcD5rI2YkISM8NyQkIiRnJiEiIiQiMjpOSnQzKWZ6TiEjPDckJCIkcSYhIiIkIjE9dztBUCRSYyQhIzs3JCQiJCFlISIiJCIxUWZYIVtVY2EkISM7NyQkIiQhZiEiIiQiMnVtUG0qb3JDTiEjPDckJCIkKychIiIkIjE4IlwqSCVbNl0kISM7NyQkIiQ1JyEiIiQiMUJjJDRZRlxaJCEjOzckJCIkPychIiIkIjI7VER2R1ZnVyQhIzw3JCQiJEknISIiJCIydF5oMTclWzlNISM8NyQkIiRTJyEiIiQiMlY1IVx6cUIhUSQhIzw3JCQiJF0nISIiJCIxYEo1Jz0pR1ZMISM7NyQkIiRnJyEiIiQiMSZRTClmQmkuTCEjOzckJCIkcSchIiIkIjEoNCo0Nk1BaEshIzs3JCQiJCFvISIiJCIyY1Inek1TMjtLISM8NyQkIiQhcCEiIiQiMUdoQi5lOm9KISM7NyQkIiQrKCEiIiQiMUxFaWUiXHU2JCEjOzckJCIkNSghIiIkIjI4Q3JdU0xSMSQhIzw3JCQiJD8oISIiJCIyJGU3OytuZTJJISM8NyQkIiRJKCEiIiQiMi8rQGcvJ1FbSCEjPDckJCIkUyghIiIkIjFENyUqenNJJylHISM7NyQkIiRdKCEiIiQiMlc4OFgxRDgjRyEjPDckJCIkZyghIiIkIjFLcEd4R1RgRiEjOzckJCIkcSghIiIkIjIkKnBZKioqSGEjbyMhIzw3JCQiJCF5ISIiJCIxJT0oKnBdJ28zRSEjOzckJCIkIXohIiIkIjJaKHkmUkQ4PWAjISM8NyQkIiQrKSEiIiQiMWQmZl0nPSo9WCMhIzs3JCQiJDUpISIiJCIxMCRRL3MqKSlvQiEjOzckJCIkPykhIiIkIjIxckBHJUh4I0cjISM8NyQkIiRJKSEiIiQiMSUqPXEkUTFOPiMhIzs3JCQiJFMpISIiJCIxaicqKSozTzAsQCEjOzckJCIkXSkhIiIkIjJKJ1I5JSlvUDA/ISM8NyQkIiRnKSEiIiQiMiVbMHhlclYxPiEjPDckJCIkcSkhIiIkIjJvZSRvXVM+Lz0hIzw3JCQiJCEpKSEiIiQiMSpwQiVIZWcpcCIhIzs3JCQiJCEqKSEiIiQiMmJwRCMqUkgnKmUiISM8NyQkIiQrKiEiIiQiMjAySDlHP3NaIiEjPDckJCIkNSohIiIkIjFqU0snZks4TyIhIzs3JCQiJD8qISIiJCIyQCllVVchPj5DIiEjPDckJCIkSSohIiIkIjJRWHByKTMkKj02ISM8NyQkIiRTKiEiIiQiMWMhbyxGekojKiohIzw3JCQiJF0qISIiJCIxWygqb3BbRz8nKSEjPDckJCIkZyohIiIkIjB0U1oqUTQhRyghIzs3JCQiJHEqISIiJCIyQy1Ic01mPyFmISM9NyQkIiQhKSohIiIkIjJtaktLeT1jWyUhIz03JCQiJCEqKiEiIiQiMjpXUU0zJT5JSSEjPTckJCIlKzUhIiIkIjJGOCQ9PDc+TjohIz0tJiUmX0FYSVNHNiMiIiI2Jy0lK19HUklETElORVNHNictJSZDT0xPUkc2JiUkUkdCRyQiIiEhIiIkIiIhISIiJCIiISEiIi0lKkxJTkVTVFlMRUc2IyIiIS0lKlRISUNLTkVTU0c2IyIiIS0lLVRSQU5TUEFSRU5DWUc2IyQiIiEhIiItJSlfVklTSUJMRUc2IyIiIS0lJkNPTE9SRzYmJSRSR0JHJCIiISEiIiQiIiEhIiIkIiIhISIiLSUqTElORVNUWUxFRzYjIiIhLSUqVEhJQ0tORVNTRzYjIiIhLSUtVFJBTlNQQVJFTkNZRzYjJCIiISEiIi0mJSZfQVhJU0c2IyIiIzYnLSUrX0dSSURMSU5FU0c2Jy0lJkNPTE9SRzYmJSRSR0JHJCIiISEiIiQiIiEhIiIkIiIhISIiLSUqTElORVNUWUxFRzYjIiIhLSUqVEhJQ0tORVNTRzYjIiIhLSUtVFJBTlNQQVJFTkNZRzYjJCIiISEiIi0lKV9WSVNJQkxFRzYjIiIhLSUmQ09MT1JHNiYlJFJHQkckIiIhISIiJCIiISEiIiQiIiEhIiItJSpMSU5FU1RZTEVHNiMiIiEtJSpUSElDS05FU1NHNiMiIiEtJS1UUkFOU1BBUkVOQ1lHNiMkIiIhISIiLSUpX1ZJU0lCTEVHNiMiIiItJSVST09URzYnLSUpQk9VTkRTX1hHNiMkIiQhUSEiIi0lKUJPVU5EU19ZRzYjJCIjcSEiIi0lLUJPVU5EU19XSURUSEc2IyQiJV1KISIiLSUuQk9VTkRTX0hFSUdIVEc2IyQiJVNMISIiLSUpQ0hJTERSRU5HNiItJStBTk5PVEFUSU9ORzYnLSUpQk9VTkRTX1hHNiMkIiIhISIiLSUpQk9VTkRTX1lHNiMkIiIhISIiLSUtQk9VTkRTX1dJRFRIRzYjJCIlK1MhIiItJS5CT1VORFNfSEVJR0hURzYjJCIlK1MhIiItJSlDSElMRFJFTkc2Ig==Ig==