> | 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]); |
> |