中国科学院大学数值分析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。