數值分析 最小二乘 matlab
1. 已知函式在下列各點的值為
| -1 |
-0.75 |
-0.5 |
0 |
0.25 |
0.5 |
0.75 |
|
| 1.00 |
0.8125 |
0.75 |
1.00 |
1.3125 |
1.75 |
2.3125 |
分別用一次、二次、三次最小二乘擬合多項式擬合上述資料，畫出所給資料和所求最小二乘擬合多項式的影象。
程式：
function f=multifit(x,y,wfunc,n)
syms t
%x,y為給定資料陣列，wfunc為權函式，n為要求擬合多項式的次數
N=length(x);
M=length(y);
if(N ~= M)
disp('x與y維數不匹配');
return;
end
var = findsym(sym(wfunc));
w = subs(wfunc,'var',x);
g(1:(2*n+1))=0;
b(1:(n+1))=0;
for j=1:(2*n+1)
for k=1:N
g(j)=g(j)+w(j)*x(k)^(j-1);
if(j<(n+2))
b(j)=b(j)+w(j)*y(k)*x(k)^(j-1);
end
end
end
G(1,:)=g(1:(n+1));
for i=2:(n+1)
G(i,:)=g(i:(n+i));
end
coff=b'\G;
f = coff(1);
l = 1;
for i=1:n
l = l*t;
f = f+coff(i+1)*l;
end
一維：
x=[-1 -0.75 -0.5 0 0.25 0.5 0.75];
y=[1.00 0.8125 0.75 1.00 1.3125 1.75 2.3125];
plot(x,y)
hold on
w=ones(1,7);
syms t
f=multifit(x,y,w,1)
t=-1:0.02:1;
tf=subs(f,'t',t);
plot(t,tf)
f =
1723890370956185/2251799813685248 - (923516587282913*t)/18014398509481984
二維：
f =
(2832971806309577*t^2)/9007199254740992 - (5864710038001133*t)/72057594037927936 + 6907316203181555/9007199254740992
三維：
f =
- (8450352688460543*t^3)/72057594037927936 + (89154935708507*t^2)/281474976710656 - (6164715222057551*t)/72057594037927936 + 6924830714825963/9007199254740992
2. 已知一組實驗資料如下，
| 2 |
3 |
6 |
9 |
10 |
|
| -0.5 |
1.2 |
3.1 |
4.5 |
7.3 |
|
| 0.125 |
0.125 |
0.25 |
0.125 |
0.375 |
求擬合上述資料的二次最小二乘擬合多項式。
主程式：
x=[2 3 6 9 10];
y=[-0.5 1.2 3.1 4.5 7.3];
w=[0.125 0.125 0.25 0.125 0.375];
plot(x,y)
hold on
syms t
f=multifit(x,y,w,2)
t=2:0.1:10;
tf=subs(f,'t',t);
plot(t,tf)
legend('已知點','擬合')
影象：
f =
(8318828653518565*t^2)/562949953421312 + (7344269211390441*t)/4503599627370496 + 6836208550152757/36028797018963968
3. 設，分別求次數為2,3,6,8的多項式，使得
達到最小，並畫出和的曲線進行比較。
程式：
function f = Legendre(func,n)
%求函式在[-1,1]的關於權函式為1的n次最佳平方逼近多項式f，並計算插值多項式f在資料點x0的函式值f0
syms t;
P(1:n+1) = t;
P(1) = 1;
P(2) = t;
c(1:n+1) = 0.0;
c(1)=int(subs(func,findsym(sym(func)),sym('t'))*P(1),t,0,1)/2;
c(2)=3*int(subs(func,findsym(sym(func)),sym('t'))*P(2),t,0,1)/2;
f = c(1)+c(2)*t;
for i=3:n+1
P(i) = ((2*i-3)*P(i-1)*t-(i-2)*P(i-2))/(i-1);
c(i) = (2*i-1)*int(subs(func,findsym(sym(func)),t)*P(i),t,0,1)/2;
f = f + c(i)*P(i);
if(i==n+1)
f = vpa(f,6);
end
end
n=2:
syms x
fun=sin(pi*x);
t=0:0.1:1;
tf=subs(fun,'x',t);
plot(t,tf)
hold on
f = Legendre(fun,2)
f2=subs(f,'t',t);
plot(t,f2)
f =
- 0.128828*t^2 + 0.477465*t + 0.361253
N=3:
f =
0.863545*t - 0.19304*t*(7.5*t^2 - 2.5) - 0.128828*t^2 + 0.361253
N=6:
f =
0.922023*t - 0.222279*t*(7.5*t^2 - 2.5) - 0.021929*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375) - 0.0383952*t*(2.93333*t*(7.5*t^2 - 2.5) - 5.86667*t + 2.2*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375)) + 0.144211*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 0.52012*t^2 + 0.144936
N=8:
f =
0.925827*t - 0.224181*t*(7.5*t^2 - 2.5) - 0.0233554*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375) - 0.0568912*t*(2.93333*t*(7.5*t^2 - 2.5) - 5.86667*t + 2.2*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375)) + 0.0158537*t*(6.85714*t - 3.42857*t*(7.5*t^2 - 2.5) - 2.57143*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375) + 2.14286*t*(2.16667*t*(2.93333*t*(7.5*t^2 - 2.5) - 5.86667*t + 2.2*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375)) - 2.70833*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) - 12.1875*t^2 + 4.0625)) + 0.167331*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 0.62416*t^2 + 0.00118873*t*(2.16667*t*(2.93333*t*(7.5*t^2 - 2.5) - 5.86667*t + 2.2*t*(2.25*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) + 10.125*t^2 - 3.375)) - 2.70833*t*(4.66667*t - 2.33333*t*(7.5*t^2 - 2.5)) - 12.1875*t^2 + 4.0625) + 0.110256