Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- //Феномен Рунге
- function y=my_func(x)
- y = 1./(1 + x.^2)
- endfunction
- function y=newton_method(xCount, x, xInterpolCount, Xn, Yn)
- //заполняем нулями и высчитываем диагональ
- for i=1:xInterpolCount
- for j=1:xInterpolCount
- if i==j then
- M(i,j) = Yn(i);
- else
- M(i,j) = 0;
- end;
- end;
- end;
- //высчитываем разделённые разности
- for i=1:+1:xInterpolCount-1
- for j=1:xInterpolCount-i
- M(j,j+i)=(M(j+1,j+i) - M(j,j+i-1)) / (Xn(j+i)-Xn(j));
- end
- end;
- for i = 1:xCount
- y(i) = Yn(1);
- if (xInterpolCount > 1)
- then
- p=1;
- for u=2:+1:xInterpolCount
- p = p*(x(i) - Xn(u-1));
- y(i) = y(i) + (M(1,u)*p);
- end
- end
- end
- endfunction
- function p=PolyLagr(m,n,x,xt,ft)
- // m кол-во Х n- кол-во узлов
- //x-наши задаваемые Х xt-узлы
- // ft - знач в узлах
- for k=1:m
- sum=0
- for i=1:n
- w(i)=1
- for j=1:n
- if (j ~= i)
- w (i) =w(i) * (x (k) - xt (j)) / (xt (i) - xt (j));
- end
- end
- sum = sum + ft(i) * w(i);
- end
- p(k) =sum;
- end
- endfunction
- function RP=RealP(n,f,f1)
- // f-данная функция f1- интерполяц. полином
- for i=1:n
- RP(i)= abs(f(i)-f1(i))
- end
- endfunction
- // M- максимум в 10 производной
- function F=FormulaP()
- endfunction
- xCount = 100;
- xInterpolCount = 10;
- a = -5;
- b = 5;
- h1 = (b-a)/(xCount - 1); //шаг для иксов функции
- h2 = (b-a)/(xInterpolCount - 1); //шаг для узлов функции
- x1 = [a:h1:b]; //иксы функции
- x2 = [a:h2:b]; //узлы интерполяции
- f1 = my_func(x1); //функция в иксах
- F1 = my_func(x2); //функция в узловых точках
- fN = newton_method(xCount, x1, xInterpolCount, x2, F1);//полином Ньютона
- fL=PolyLagr(xCount,xInterpolCount,x1,x2,F1)
- RealN=RealP(xCount,f1,fN)
- RealL=RealP(xCount,f1,fL)
- subplot(3,4,1)
- plot(x1,f1,"r",x1,fN,"b");
- xtitle('кр.-исходная функция, Син.-Ньютон')
- subplot(3,4,2)
- plot(x1,f1,"r",x1,fL,"g");
- xtitle('кр.-исходная функция, зел.-Лагранж')
- subplot(3,4,3)
- plot(x1,RealN,"g");
- xtitle('реал. погр. пол.Ньютона')
- subplot(3,4,4)
- plot(x1,RealL,'r');
- xtitle('реал. погр. пол.Лагранжа')
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement