Блок схема к программе

> restart:
> Digits:=25:
> with(linalg):
> N:=7:
> l:=1:
> lambda:=(i)->(evalf(Pi)*i/l)^2:        собственные числа оператора Лапласа
> v:=(i,x)->sin(evalf(Pi)*i/l*x):     собственные функции оператора Лапласа
> 
> q:=(x)->exp(2*x+1):                 потенциал
> A:=matrix(N,N):                     матрица оператора Штурма-Лиувилля
> for i from 1 to N do
>   for j from 1 to N do
>     A[i,j]:=int((diff(v(i,x),x,x)+q(x)*v(i,x))*v(j,x),x=0..l);
>   od:
> od:
> Y0:=matrix(N,1,0):          нулевое приближение
> Y0[1,1]:=1:  
> LAMBDA:=array(1..N):        собственные числа оператора Штурма-Лиувилля
> X:=matrix(N,N):             собственные функции оператора Штурма-Лиувилля 
> Y:=matrix(N,N):             основная матрица СЛАУ для поиска коэффициентов хар-го ур-ния
> B:=array(1..N):             столбец свободных членов СЛАУ для поиска коэффициентов хар-го ур-ния
> P:=array(1..N):             коэф-ты хар-го уравнения
> Q:=array(0..N-1,1..N):      коэф-ты произвольно многочлена, вычисляемы по схеме Горнера
>  
> C:=matrix(N,N):             вспомогательная матрица
> nevyazka:=array(1..N):      массив значений невязки полученного решения
> 
     составляем и решаем систему уравнений для поиска коэф-тов хар-го полинома
> for j from 1 to N do
>   for i from 1 to N do 
>     Y[i,N-j+1]:=Y0[i,1]:
>   od:
>   Y0:=evalm(A&*Y0): 
> od:
> B:=evalm(-Y0):
> P:=linsolve(Y,B):
> 
       составляем характирестический полином и находим его корни
>  k:=lambda^N:                
>  for i from 1 to N do
>    k:=k+P[i,1]*lambda^(N-i):
>  od:
>  LAMBDA:=solve(k=0,lambda):  
> 
        находим коэф-ты произвольного многочлена по схеме Горнера
>  for i from 1 to N do         
>    Q[0,i]:=1:
>    for j from 1 to N-1 do 
>      Q[j,i]:=evalf(LAMBDA[i]*Q[j-1,i]+P[j,1]):
>    od:
>  od:
> 
       находим собственные функции
>  for k from 1 to N do         
>    for i from 1 to N do
>      X[i,k]:=Y[i,1]:
>      for j from 2 to N do
>        C[i,j]:=Y[i,j]*Q[j-1,k]:
>        X[i,k]:=X[i,k]+C[i,j]
>      od:
>    od:
>  od:
> 
       нормируем собственные функции
>  for j from 1 to N do        
>    delta:=0:
>    for i from 1 to N do
>      delta:=delta+X[i,j]^2:
>    od:
>    delta:=delta^(1/2): 
>    for i from 1 to N do
>      X[i,j]:=X[i,j]/delta:
>    od:
>   выполним проверку, посчитав невязку полученного решения
>   E:=matrix(N,N,0):          
>   for i from 1 to N do
>     E[i,i]:=LAMBDA[j]:
>   od:
>   X1:=matrix(N,1):
>   for i from 1 to N do
>     X1[i,1]:=X[i,j]:
>   od:
>   K:=evalf(evalm((A-E)&*X1)):
>   nevyazka[j]:=0:
>   for i from 1 to N do
>     nevyazka[j]:=nevyazka[j]+K[i,1]^2:
>   od:
>   nevyazka[j]:=nevyazka[j]^(1/2):
> od:  
> 
> for i from 1 to N do
>    lprint(`    i=`,i,`     LAMBDA[i]=`,LAMBDA[i],`    nevyazka=`,nevyazka[i]); 
> od: