> N3:=x->seq(subs(x0=a+(i-1)*h,x1=a+(i-1)*h+h/2, x2=a+i*h,(x-x0)*(x-x1)/((x2-x0)*(x2-x1))),i=1..n);
Есть другой, более универсальный, способ задания базисных функций. Он основан на вычислении коэффициентов {Cie }
в выражении типа (10.5) из каких-либо условий, например, что функция в узлах {xi} должна принимать значения { i}. Кстати, эти условия могут быть наложены и на производные, и здесь мы приходим к эрмитовой интерполяции, сплайн-интерполяции и др. Разберите следующий листинг программы, реализующий ввод этим универсальным способом базисных функций интерполяции Лагранжа 2-го порядка. (Кстати, этот способ лучше приспособлен для случая неравноотстоящих узлов.)
>restart;
>a:=1;b:=2;n:=4;h:=(b-a)/n; m:=2*n;
>wt:=x->c1*x^2+c2*x+c3;
>wr:=solve({wt(xi)=w1,wt(xi+h/2)=w2,wt(xi+h)=w3},
{c1,c2,c3});
>uf:=collect(subs(wr,wt(x)),{w1,w2,w3});
uv:=[diff(uf,w1),diff(uf,w2),diff(uf,w3)];
>N1:=x->seq(subs(xi=a+(i-1)*h,uv[1]),i=1..n): N2:=x->seq(subs(xi=a+(i-1)*h,uv[2]),i=1..n): N3:=x->seq(subs(xi=a+(i-1)*h,uv[3]),i=1..n):
Следующий шаг – задание аппроксимации неизвестной функции. Задавая узловые величины через массив y, имеем
> z:=x->seq(y[2*i-2]*N1(x)[i]+y[2*i-1]*N2(x)[i]+ y[2*i]*N3(x)[i],i=1..n);
Дальнейшие шаги почти полностью повторяют строки [4]–[8] листинга из пункта 7. Следует только учесть, что переменная m – число узлов минус один – здесь определяется иначе.
> F:=0; for i from 1 to n do
F:=F+(1/2)*int(x^2*diff(z(x)[i],x)^2+ 2*expand((z(x)[i]^2))-2*(1+2./x)*z(x)[i],
x=a+(i-1)*h..a+i*h);od: |
# |
[4] |
> F:=F-4*z(2)[n]: |
# |
[5] |