ЧМ ЛР3
.docx
|
t |
0 |
1/4 |
1/2 |
3/4 |
1 |
|
f(t) |
0 |
0.7071 |
0.5 |
0.7071 |
0 |
>> X = [0, 1/4, 1/2, 3/4, 1];
>> X1 = [0, 1/6, 1/3, 1/2, 1];
plot(X,sin(pi.*X),'-b');
grid on
hold on
S = 0:pi/100:1;
plot(S, sin(pi.*S), '-k')
Y(1) = sin(pi*X1(1));
for i = 2:1:length(X)
Y(i) = sin(pi*X(i-1)) + ((sin(pi*X(i))-sin(pi*X(i-1)))/(X(i)-X(i-1)))*(X1(i)-X(i-1));
end
plot(X1,Y,'or')
Y = 0 0.4714 0.8047 1.0000 0.0000
>> sin(pi.*X1)
ans = 0 0.5000 0.8660 1.0000 0.0000
function Lagr(X, X1, f)
syms x
for i=1:1:length(X)
F(i)=subs(f,X(i));
end
P = 0;
p = 0;
for i = 1:1:length(X)
A = 1;
for j = 1:1:length(X)
if (j ~= i)
A = A * ((x - X(j))/(X(i) - X(j)));
end
end
p = p + F(i)*A;
end
P = subs(p, X1);
plot(X1, P, '-k')
vpa(P,4)
Y = 0:pi/100:1;
R = abs(subs(f, Y) - subs(p, Y));
A = max(R);
vpa(A)
end
>> X = [0, 1/4, 1/2, 3/4, 1];
>> X1 = [0, 1/6, 1/3, 1/2, 1];
>> ezplot('sin(pi*x)',[0,1])
>> hold on
>> grid on
>> Lagr(X,X1,@(x)sin(pi*x))
ans =[ 0, 0.4991, 0.8663, 1.0, 0]
По линейной интерполяции: Y = 0 0.4714 0.8047 1.0000 0.0000
Реальное значение: ans = 0 0.5000 0.8660 1.0000 0.0000
ans = 0.0017873298346041081364927178839716
function Newton(X, X1, f)
for i=1:1:length(X)
F(i)=subs(f,X(i));
end
A(1) = F(1);
A(2) = (F(2) - A(1))/(X(2)-X(1));
for i=3:1:length(F)
B = A(1);
for j = 2:1:i-1
C = 1;
for k = 1:1:j-1
C = C * (X(i) - X(k));
end
B = B + A(j)*C;
end
A(i) = (F(i) - B)/(C * (X(i) - X(i-1)));
end
syms x
p = A(1);
P(1) = A(1);
for i = 2:1:length(X)
C = 1;
for j = 1:1:i-1
C = C * (x - X(j));
end
p = p + A(i)*C;
end
P = subs(p, X1);
plot(X1, P)
P = vpa(P,4)
R = F - P;
R = vpa(R,4)
end
>> X = [0, 1/4, 1/2, 3/4, 1];
>> X1 = [0, 1/6, 1/3, 1/2, 1];
>> ezplot('sin(pi*x)',[0,1])
>> hold on
>> grid on
>> Newton(X,X1,@(x)sin(pi*x))
P = [ 0, 0.4991, 0.8663, 1.0, 0]
По линейной интерполяции: Y = 0 0.4714 0.8047 1.0000 0.0000
Реальное значение: ans = 0 0.5000 0.8660 1.0000 0.0000
По Лагранжу: ans = [ 0, 0.4991, 0.8663, 1.0, 0]
ans = 0.0017873298346041081364927178839716
Смотреть функции Lagr() и Newton() в предыдущем пункте
function Lagr(X, X1, f)
for i=1:1:length(X)
F(i)=subs(f,X(i));
end
for k = 1:1:length(X1)
A = ones(length(X1), length(X1));
for i = 1:1:length(X)
for j = 1:1:length(X)
if (j ~= i)
A(i,j) = ((X1(k) - X(j))/(X(i) - X(j)));
end
end
end
P(k) = F*prod(A,2);
end
plot(X1, P, '-k')
vpa(P, 4)
end
>> X = [0, 1/4, 1/2, 3/4, 1, 2.1];
>> X1 = [0, 1/6, 1/3, 1/2, 1, 2];
>> Lagr(X,X1,@(x)sin(pi*x))
ans = [ 0, 0.4961, 0.868, 1.0, 0, -0.1286]
>> X = -5:1:5;
>> f = @(x)1/(1+25*x^2);
>> X1 = [-4.5, -3.5, -2.5, -1.5, -0.5, 0, 0.5, 1.5, 2.5, 3.5, 4.5];
>> Lagr(X,X1,f)
ans = [ 4.65, -0.9391, 0.3932, -0.28, 0.68, 1.0, 0.68, -0.28, 0.3932, -0.9391, 4.65]
>> grid on
>> hold on
>> ezplot('1/(1+25*x^2)',[-5,5])
>> axis equal
