华中科技大学数值分析-数值分析实验程序

华中科技大学数值分析-数值分析实验程序

Function t_charpt1_1

%数值试验1.1病态问题

%输入：[0 20]之间的扰动项及小的扰动常数

%输出：加扰动后得到的全部根

Result=inputdlg({‘请输入扰动项：在[0 20]之间的整数:’},’charpt1_1’,1’{‘19’});

Numb=str2num(char(result));

if ((Numb>20)|(Numb<0))

errordlg(‘请输入正确的扰动项：[0 20]之间的整数!’);

return;

end

result=inputdlg({‘请输入(0 1)之间的扰动常数：’},’charpt1_1’,1,{‘0.00001’}); ess=str2num(char(result));

ve=zeros(1,21);

ve(21-Numb)=ess;

root=roots(poly(1:20)+ve);

dis p([‘对扰动项’,num2str(Numb),’加扰动’,num2str(ess),’得到的全部根为:’]); disp(num2str(root));

function charpt3

%数值实验三：含“实验3.1”和“实验3.2”

%子函数调用:dlsa

%输入：实验选择

%输出：原函数及求得的相应插值多项式的函数的图像以及参数alph和误差r

result=inputdlg({‘请选择实验，若选 3.1，请输入1，否则输入2:’},’charpt_3’,1,{‘1’});

Nb=str2num(char(result));

if(Nb~=1&(Nb~=2)errordlg(‘实验选择错误！’))；

return;

end

x0=-1:0.5:2

y0=[-4.447 -0.452 0.551 0.048 -0.447 0.549 4.552];

n=3;%n为拟合阶次

if(Nb==1)

alph=polyfit(x0,y0,n);

y=polyval(alph,x0);

r=(y0-y)*(y0-y)’;%平方误差

x=-1:0.01:2;

y=polyval(alph,x);

plot(x,y,’k-’);

xlabel(‘x’);ylabel(‘y0*and

ployfit.y-’);

hold on;

plot(x,y,’k-’);

title(‘离散数据的多项式拟合’);

grid on;

else

result=inpurdlg({‘请输入权向量w:’},’charpt_3,1,{‘[1 1 1 1 1 1 1]’});

w=str2num(char(result));

[a,b,c,alph,r]=dlsa(x0,y0,w,n);

end

disp([‘平方误差:’,sprint(‘%g,r’)]); disp([‘参数a lph:’,sprint(‘%\t’,alph)])

%---------------------------------------------------------------------------------

---------------------------------------------function[a,b,c,alph,r]=dlsa(x,y,w,n) %功能：用正交化方法对离散数据作多项式最小二乘拟合。

%输入：m+1个离散点（x,y,w）,x,y,w分别用行向量给出。

% 拟合多项式的次数n,0

%输出:三项递推公式的参数a,b,拟合多项式s(x)的系数c和alph，

% 平方误差r=(y-s,y-s),并作离散点列和拟合曲线的图形

m=length(x)-1;

if(n<1|n>=m)

errordlg(‘错误：n<1或者n>=m!’); return;

end

%求三项递推公式的参数a,b,拟合多项式s(x)的系数c,其中d(k)=(y,sk);

s1=0; s2=ones(1,m+1); v2=sum(w);

d(1)=y*w’;c(1)=d(1)/v2;

for k=1:n

xs=x.*s2.^2*w’; a(k)=xs/v2;

if(k==1)

b(k)=0;

else

b(k)=v2/v1;

end

s3=(x-a(k)).*s2-b(k)*s1;

v3=s3.^2*w’;

d(k+1)=y.*s3*w’;c(k+1)=d(k+1)/v3 end

%求平方误差r

r=y.*y*w’-c*d’;

%,求拟合多项式s(x)的降幂系数alph

alph=zeros(1,n+1);T=zeros(n+1,n+2);

T(:,2)=ones(n+1,1);T(2,3)=-a(1);

if(n>=2)

for k=3:n+1

for i=3:k+1

T(k,i)=T(k-1,i)-a(k-1)*T(k-1,i-1)-b(k-1)* T(k-2,i-2);

end

end

end

for i=1:n+1

for k=i:n+1

alph(n+2-i)=alph(n+2-i)+c(k)*T(k,k+2-i); end

end

%用秦九韶方法计算s(t)的输出序列(t,s)

xmin=min(x); xmax=max(x); dx=(xmax-xmin)/(25*m);

t=(xmin-dx):dx:(xmax+dx);

s=alph(1);

for k=2:n+1

s=s.*t+alph(k);

end

%输出点列x-y和拟合曲线t-s的图形

plot(x,y,’*’,t,s,’-’);

title(‘离散数据的多项式拟合’)；

x label(‘x’);ylbel(‘y’);

grid on;

function charpt5_1

%数值实验5.1：常微分方程性态和R-K法稳定性实验

%输入：参数a，步长h

%输出：精确解和数值解图形对比

clf;

result=inputdlg({‘请输入[-50 50]间的参数a:’},’实验5.1’,1,{‘40’});

a=str2num(char(result));

if((a<-50)|(a>50))

errordlg(‘请输入正确的参数a！’)；

return;

end

result=inputdlg({‘请输入（0 1）之间的步长：’}，‘实验5.1’，1，{‘0.01’})；h=str2num(char(result));

if((h>=1)|(h<=0))

errordlg(‘请输入正确的（0 1）间的步长！’)；

return；

end

x=0:h:1;

y=x;

N=length(x);

y(1)=1;

func=inline(‘1+(y-x).*a’);

for n=1:N-1

k1=func(a,x(n),y(n));

k2=func(a,x(n)+h/2,y(n)+k1*h/2); k3=func(a,x(n)+h/2,y(n)+k2*h/2); k4=func(a,x(n)+h/2,y(n)+k3*h/2); y(n+1)=y(n+h*(k1+2*k2+2*k3+k4)/6; end

y0=exp(a*x)+x;

plot(x,y0,’g+’);

hold on;

plot(x,y,’b--’);

xlabel(‘x’);

ylabel(‘y0=exp(ax)+x:+and R-K(x)--’)