我在做一个数值分析的题目中遇到一个这样的问题。题目如下:
上机用正交多项式拟合曲线sin(x),要求:y=sin(x),x在0~pi/2,取步长h=1度,做一张函数表。权函数w恒等于1。
1、输出:[S(xi)-yi]^2的总和
2、若要使精确度为10^-5,问要做多少次拟合。
由于这道题是老师给我们出的课外题目,并没有标准答案和现成的程序。我后来决定用Matlab来做,因为Matlab里面有符号工具箱,而要求的构造的正交多项式和要求的拟合多项式都是函数。我已经把程序实现了,并且画出拟合图像,可以看出4次的拟合效果非常好。我把程序附上。(为了方便大家阅读程序,我在程序中做了注释。)下面我提出我的问题:
1、我把所求的拟和函数S(x)和正交函数族P[1]~P[5]求出来了,但是输出的时候格式很不好,那些多项式的系数和常数项一输出来和就是一个长长的分数(分子和分母都是很大的数),我想让它们都变成位数可以控制的分数,请问应该这么做。我用pretty(),simplify()结果更加恶心。。。
2、由于求出的结果是很长很长的一串式子(其实不长,如果能把分数化成小数的话就是一个四次多项式,并不长。),所以在画图的时候也很耗时,4次拟合还能画,如果是次数更高就画不了,直接报错了。所以第一个问题一定得解决。
3、由于我用了符号工具箱做计算,所以程序运行起来很耗时,我输入拟合次数为4的时候,运行了将近18秒(我的机子是1.9G CPU,256内存),如果输入的次数更高就让人等得不耐烦了。请帮帮忙提出更好的决绝方案。
我迫切等待各位的热心解答。谢谢:)
-------------------------------------------------------------------------------
源程序(zhengjiao2.m)如下:(在本帖子后我把源程序作为附加上传了,下载后把txt后缀改为.m就可以在Matlab下运行。)
function zhengjiao2()
%最小二乘法曲线拟合[0,pi/2]上的sin(x)
%先生成0到90度均匀的91个节点,然后构造正交多项式做最小二乘拟合
%2006-10-29 0:21
clear;
n=input('Input the times of poly:n=');
tic;
X=0:pi/180:pi/2;
m=length(X);
Y=zeros(m,4); %Y矩阵有四列,第一列是用来存储拟合值;第二列储存真实值;
%第三列存储拟合值与真实值的差的绝对值;第四列存储差值的二范数。
Y(:,2)=sin(X);
syms x P sum; %P正交多项式族,sum是要求的拟合曲线
P(1)=1;
P(n+1)=0; %把P数组的长度定好
sum=0;
sumwucha=0; %Y第四列的总和,即是误差二范数的总和
for k=1:n
pai1=0;
tai1=0;
pai2=0;
tai2=0; %我有一个习惯,分子用pai表示,分母用tai表示
if k==1
for i=1:m
pai1=pai1+X(i)*subs(P(k),X(i))^2;
tai1=tai1+subs(P(k),X(i))^2;
end
P(k+1)=(x-pai1/tai1)*P(k);
else
for i=1:m
pai1=pai1+X(i)*subs(P(k),X(i))^2;
tai1=tai1+subs(P(k),X(i))^2;
pai2=pai2+subs(P(k),X(i))^2;
tai2=tai2+subs(P(k-1),X(i))^2;
end
P(k+1)=(x-pai1/tai1)*P(k)-pai2/tai2*P(k-1);
end
end %构造正交多项式族
for k=1:n+1
pai3=0;
tai3=0;
for i=1:m
pai3=pai3+Y(i,2)*subs(P(k),X(i));
tai3=tai3+subs(P(k),X(i))^2;
end %求出各个不同次数的正交多项式的系数
sum=sum+pai3/tai3*P(k);
end %求出拟合的多项式S(x)
sprintf('S(x)=%s',char(sum))
ezplot(sum,[0,pi/2])
grid on;
for i=1:m
Y(i,1)=subs(sum,X(i)); %各个点的拟合值
Y(i,3)=abs(Y(i,1)-Y(i,2)); %计算各个点处的误差
Y(i,4)=Y(i,3)^2; %计算各个点处误差的2范数
sumwucha=sumwucha+Y(i,4); %计算误差2范数的总和
end
Y
sumwucha
for i=1:n+1
sprintf('P[%d]=%s',i,char(P(i)))
end %输出所构造出来的所有正交多项式族
toc
[此贴子已经被作者于2006-11-2 20:21:06编辑过]