中国科学院大学数值分析MATLAB大作业(计算机学院开设)

7.编程实现题

%{

在微电机设计计算中需要查磁化曲线表，通常给出的表是磁密B每间隔100高斯磁路每厘米长所需安匝数at的值，下面要解决B从4000至11000区间的查表问题。为节省计算机存储单元，采用每500高斯存入一个at值，在利用差分公式计算。从差分表中看到三阶差分近似于0，计算时只需两阶差分。当4000≤B ≤10500时用牛顿前插公式；当10500≤B≤11000时用牛顿后插公式；试在计算机上编程实现求任一在区间［4000，11000］内的函数值。

%}

syms y0

y0=[1.38,1.48,1.58,1.69,1.81,1.94,2.10,2.28,2.50,2.76,3.06,3.41,3.83,4.33,4.93];

b=input('enter the number of b:');

if b>10500&&b<=11000

t=(b-11000)/500.0;

r=4.93+0.6*t+0.1*t*(t+1)/2;

disp(r)

elseif b<=10500&&b>=400

k=(b-4000)/500;

m=floor(k);

t=k-m;

r=y0(m+1)+(y0(m+2)-y0(m+1))*t+t*(t-1)*(y0(m+3)-2*y0(m+2)+y0(m+1))/2;

disp(r)

else

disp('the number is out of consideration')

end

clear

15.用正交多项式(格拉姆－施密特)作最小二乘拟合的程序

syms alpha;

syms beta;

syms a;

syms p;

x=input('enter the value of x(for example:[1,2,3,4,5]):');

y=input('enter the value of y(for example:[1,2,3,4,5]):');

w=input('enter the value of weight(for example:[0.1,0.2,0.3,0.4]):');

n=input('enter the value of n:');

m=length(x);

alpha(2)=sum(w*x)/sum(w);

a(1)=sum(w*y)/sum(w);

p2x= x-alpha(2);

alpha(3)=sum(w*x*(subs(p2x,x))^2)/sum(x*( subs(p2x,x))^2);

beta(2)=sum(w*( subs(p2x,x))^2)/sum(w);

a(2)=sum(w*y*( subs(p2x,x)))/sum(w*( subs(p2x,x))^2);

for k=3:m+1

alpha(k+1)=sum(w*x*((subs(p,{k,x},[k,x]))^2)/sum(w*( subs(p,{k,x},[k,x]))^2);

beta(k)=sum(w*( subs(p,{k,x},[k,x]))^2)/sum(w*( subs(p,{k-1,x},[k-1,x]))^2);

a(k)=sum(w*y* subs(p,{k,x},[k,x]))/sum(w*( subs(p,{k,x},[k,x]))^2);

p=(x-alpha(k+1))*subs(p,{k,x},[k,x])-beta(k)*subs(p,{k-1,x},[k-1,x]);

end,

for i=1:n+1

F=0;

F=F+a(i)* subs(p,{i,t},[i,t]);

end