数值分析作业

实验一  圆周率π值的计算

方案一：
S0=sin(pi/2)^2;
P0=2*sqrt(S0);
S=S0/2*(1+sqrt(1-S0));
P=4*sqrt(S);
eps=1.0e-14; % 设置精度
n=2;
format long
while(norm(S0-S)>=eps && norm(P0-P)>=eps)
    S0=S;
    S=S0/(2*(1+sqrt(1-S0)));
    P=2^(n+1)*sqrt(S);
    n=n+1;
    fprintf('S为')
    disp(S)
    fprintf('P为')
    disp(P)
end
fprintf('π的逼近值为')
disp(P)

方案二：
S2=1;
S=S2/(sqrt(2*(1+sqrt(1-S2^2))));
P=4*S;
P2=2*S2
eps=1.0e-14; % 设置精度
n=4;
format long
while(norm(S2-S)>=eps && norm(P2-P)>=eps)
    S2=S;
    S=S2/(sqrt(2*(1+sqrt(1-S2^2))));
    P=2^(n-1)*S;
    n=n+1;
    fprintf('S为')
    disp(S)
    fprintf('P为')
    disp(P)
end
fprintf('π的逼近值为')
disp(P)
 
方案三：
a=0;b=1;c=1/sqrt(2);d=1/4;e=1;
for k=1:20
a=b;
b=(b+c)/2;
c=sqrt(c*a);
d=d-e*((b-a)^2);
e=2*e;
k=k+1;
end
f=b^2/d
g=(b+c)^2/(4*d)

实验二 牛顿(Newton)迭代方法计算电阻RB实验

Ra=8670;C=1.0e-7;T2=1.4e-4;
T1=Ra*C*log(2)
D=T1/(T1+T2)*100
f=1/(T1+T2)
 
syms f Rb ;
x0=7740; % 初值 迭代区间确定为[7740,7745]
eps=1.0e-5;
f(Rb) = T2+(Ra*Rb*C/(Ra+Rb))*log(abs((Ra-2*Rb)/(2*Ra-Rb)));
F=diff(f);
x1=x0-f(x0)/F(x0);
k=1;
while norm(x1-x0)>=eps
    x0=x1;
    x1=x0-f(x0)/F(x0);
    k=k+1;
end
x=x1 ;
fprintf('Rb为')
disp(eval(x))

实验三 超松弛(SOR)迭代方法解大规模线性方程组

n = 10; % n = 100
% D、U、L
v = ones(1,n-1);
D = -4*eye(n); 
U = diag(v,1);
L = diag(v,-1);
d=-2*ones(n-2,1);
 
b = [-3;d;-3];
eps=1.0e-6;
 
for w=1.1:0.1:1.9
    k = 0;
    M = (D+w.*L)\((1-w).*D-w.*U);
    beta = w.*((D+w.*L)\b);
   
    x0=ones(1,n)*0;
    x=M.*x0+beta;
    while(norm(x0-x)>=eps)
     x0=x;
     x=M*x0+beta;
     k=k+1;
    end
    disp("w为"+w+"时，"+"迭代次数为："+k)
    fprintf('\n---------------------------------\n')
end

 

实验四 样条插值的应用

x=[0 70 130 210 337 578 776 1012 1142 1462 1841];
y=[0 57 78 103 135 182 214 244 256 272 275];
h=eye(1,10);
a=eye(1,9);
u=eye(1,9);
d=eye(1,9);
 
% 三弯矩法
for i= 1:10
    h(i) = x(i+1)-x(i); % 算步长
end
for j= 1:9
    a(j)=h(j+1)/(h(j)+h(j+1));  
    u(j)=1-a(j);
    d(j)=6*((y(j+2)-y(j+1))/h(j+1)- (y(j+1)-y(j))/h(j))/(h(j)+h(j+1));
    
end
 
% 自然边界条件，II型样条插值
A=a(1:8);
U=u(2:9);
Y=diag(A,1)+diag(U,-1)+2*eye(9);
M=Y\d';
M1=[0;M;0];
S=eye(1,36);


for t = 1:10
 f=@(k) M1(t)*(x(t+1)-k)^3/(6*h(t)) +  M1(t+1)*(k-x(t))^3/(6*h(t)) + (y(t)-M1(t)*h(t)^2/6)*(x(t+1)-k)/h(t) + (y(t+1)-M1(t+1)*h(t)^2/6)*(k-x(t))/h(t);
 fplot(f,[x(t),x(t+1)])
 hold on
 
 for s = 1:36
     k = 50*s;
     if x(t)<= k && k <= x(t+1)
          S(s) = f(k);
     end
 end
end
 
% I型样条插值
d0=6*((y(2)-y(1))/h(1)- 1)/h(1);
d10=6*(-(y(11)-y(10))/h(10))/h(10);
D0=[d0;d';d10];
A0=[1,a];
U0=[u,1];
Y0=diag(A0,1)+diag(U0,-1)+2*eye(11);
M0=Y0\D0;
S0=eye(1,36);
 
for t = 1:10
 g=@(k) M0(t)*(x(t+1)-k)^3/(6*h(t)) +  M0(t+1)*(k-x(t))^3/(6*h(t)) + (y(t)-M0(t)*h(t)^2/6)*(x(t+1)-k)/h(t) + (y(t+1)-M0(t+1)*h(t)^2/6)*(k-x(t))/h(t);
 fplot(g,[x(t),x(t+1)])
 hold on
 
 for s = 1:36
     k = 50*s;
     if x(t)<= k && k <= x(t+1)
          S0(s) = g(k);
     end
 end
end