数值计算作业(maltlab实现)

function x  = nagauss( a,b,flag )
if nargin<3
    flag=0;
end
 n=length(b);
 a=[a,b];
for k=1:(n-1)
    a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))-a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));
    a((k+1):n,k)=zeros(n-k,1);
    if flag==0
        a
    end
end
x=zeros(n,1);
 x(n)=a(n,n+1)/a(n,n);
 for k=(n-1):-1:1
     x(k)=(a(k,n+1)-a(k,(k+1):n)*x((k+1):n))/a(k,k);
 end

A=[2 3 1;1 2 1;2 1 1]
b=[450000;300000;350000]
x  = nagauss(A,b,1)

A =

     2     3     1
     1     2     1
     2     1     1


b =

      450000
      300000
      350000


x =

      100000
       50000
      100000

function Y=lagrange(x1,y1,xx)
syms x
n=length(x1);
for i=1:n
t=x1;t(i)=[];L(i)=prod((x-t)./(x1(i)-t));
end
u=sum(L.*y1);
p = simplify(u) 
Y=subs(p,x,xx); 
end

x1= [1998 2000 2002 2004 2006 2008]
y1=[2.32 2.98 2.82 3.56 4.87 5.66]
xx = [1999 2001 2003 2005 2007]
Y=lagrange(x1,y1,xx)

x1 =

        1998        2000        2002        2004        2006        2008


y1 =

    2.3200    2.9800    2.8200    3.5600    4.8700    5.6600


xx =

        1999        2001        2003        2005        2007

 
p =
 
(43*x^5)/128000 - (64667*x^4)/19200 + (259337227*x^3)/19200 - (5200153633*x^2)/192 + (20365575359273*x)/750 - 544483909361251/50
 

Y =

2.9805    2.8242    3.0586    4.2227    5.3945

 

 

 

源文件名称	具体内容
nalupad	LU分解紧凑格式，L为单位下三角阵，U为单位上三角阵
nalu	LU分解格式，L为单位下三角阵，U为单位上三角阵
GaussJordan	高斯-约当分量形式
nagauss	顺序高斯消去法分量形式
fjacobi	Jacobi迭代分量形式
jacobi	jacobi迭代矩阵形式（未实现）
jac	jacobi迭代矩阵形式
nags	高斯-塞德尔迭代矩阵形式
gau_seidel	高斯-塞德尔迭代分量形式
nabisect	二分法解非线性（函数句柄）
sfen	二分法feval函数调用函数句柄
f	创建任意的f函数
Newton_f	牛顿迭代法
secant	双点割法
lagrange 拉格朗日插值

%拉格朗日插值
function Y=lagrange(x1,y1,xx)
syms x
n=length(x1);
for i=1:n
t=x1;t(i)=[];L(i)=prod((x-t)./(x1(i)-t));
end
u=sum(L.*y1);
p = simplify(u) 
Y=subs(p,x,xx); 
end

x1= [1998 2000 2002 2004 2006 2008]
y1=[2.32 2.98 2.82 3.56 4.87 5.66]
xx = [1999 2001 2003 2005 2007]
Y=lagrange(x1,y1,xx)


function x  = nagauss( a,b,flag )
%用途：顺序高斯消去法解线性方程组ax=b
% a:系数矩阵，b:右端列变量
% flag:若为0，则显示中间过程，否则不显示，默认值为0
% x:解向量
if nargin<3
    flag=0;
end
n=length(b);
a=[a,b];
% 消元
for k=1:(n-1)
    a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))-a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));
    a((k+1):n,k)=zeros(n-k,1);
    if flag==0
        a
    end
end
%  回代
x=zeros(n,1);
x(n)=a(n,n+1)/a(n,n);
for k=(n-1):-1:1
    x(k)=(a(k,n+1)-a(k,(k+1):n)*x((k+1):n))/a(k,k);
end

end


%列主元高斯消去
function x = nagauss(A, B)
n=length(B);
x=zeros(n,1);
c=zeros(1,n);
d1=0;
for i=1:n-1
    max=abs(A(i,i));
    m=i;
    for j=i+1:n
        if max<abs(A(j,j))
            max=abs(A(j,j));
            m=j;
        end
    end
    if(m~=i)
        for k=i:n
            c(k)=A(i,k);
            A(i,k)=A(m,k);
            A(m,k)=c(k);
        end
        d1=B(i);
        B(i)=B(m);
        B(m)=d1;
    end
    for k=i+1:n
        for j=i+1:n
            A(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i);
        end
        B(k)=B(k)-B(i)* A(k,i)/A(i,i);
        A(k,i)=0;
    end
end
x(n)=B(n)/A(n,n);
for i=n-1:-1:1
    sum=0;
    for j=i+1:n
        sum=sum+A(i,j)*x(j);
    end
    x(i)=(B(i)-sum)/A(i,i);
end

 A = [1, 2, 3; 2, 4, 9; 3, 4, 9]
 B = [4; 2; 5];
 x = nagauss(A, B)

function x=gau_seidel(A,b,P,delta,n)
%A为n维非奇异阵；b为n维值向量
%P为初值；delta为误差界；n为给定的迭代最高次数
N=length(b);
for k=1:n       
   for j=1:N             
     if j==1                 
        x(1)=(b(1)-A(1,2:N)*P(2:N))/A(1,1);           
     elseif j==N                  
         x(N)=(b(N)-A(N,1:N-1)*(x(1:N-1))')/A(N,N);
        else                  
          x(j)=(b(j)-A(j,1:j-1)*x(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);                      
     end      
   end      
   err=abs(norm(x'-P));       
   P=x'      
   if(err<delta)            
      break;       
   end
end
x=x';
k
err
end

function solution=GaussJordan(gauss,precision)
%%gauss 为用户输入的增广矩阵
%%precision 为用户输入的精度要求,如不输入或输入有误,则默认为10位
% %输入参数检查             
if nargin==2               
  try                         
   digits(precision);               
  catch                        
   disp('您输入的精度有误,这里按照缺省的精度(10位有效数字)计算')                            
   digits(10);                 
  end               
else                       
  digits(10);               
end               
A=vpa(gauss);               
row=size(A,1);               
col=size(A,2);               
if ndims(A)~=2||col-row~=1                        
   disp('矩阵大小有误,不能使用GaussJordan消去法')                    
   return              
end              
if det(gauss(:,1:row))==0                       
   disp('该方程的系数矩阵行列式为零,无解或有无穷多解,不能使用GaussJordan消去法')                        
return             
end
% %消元过程                    
for i=1:row
% %选取主元
Max=0.0;                         
    for j=i:row                               
        if double(abs(A(j,i))-Max)>0                               
           Max=abs(A(j,i));                              
           max_row=j;                       
        end                
    end                
    temp=A(i,:);                
    A(i,:)=A(max_row,:);              
    A(max_row,:)=temp;
% %上下消元               
for k=1:row                       
        if k~=i                               
           A(k,:)=vpa(A(k,:)-A(i,:)*A(k,i)/A(i,i));                        
        end                
    end       
end
% %得到解
for i=1:row          
    solution(i)=vpa(A(i,col)/A(i,i));    
end
end

function x=jac(A,b,x0,e,N)
%用途:用向量形式(普通储存格式)的jacobi迭代解线性方程组Ax=b
%格式:x=jac(A,b,x0,e,N)。A为系数矩阵，b为右端向量，x返回解向量，
%x0为初始向量(默认原点),e为精度(默认le-4),设置迭代次数上限N以防
%发散(默认500)
n=length(b);
if nargin<5,N=500;end
if nargin<4,e=1e-4;end
if nargin<3,x0=zeros(n,1);end
x=x0;x0=x+2*e;
k=0;Al=tril(A);iAl=inv(Al);
while norm(x0-x,inf)>e&&k<N,      
k=k+1;      
x0=x;
x=-inv(diag(diag(A)))*(A-diag(diag(A)))*x0+inv(diag(diag(A)))*b;      
x'
end
if k==N,warning('已达迭代次数上限');
end
end

function [ x,iter ]=jacobi( A,b,tol )
D=diag(diag(A));
L=D-tril(A);
U=D-tril(A);
x=zeros(size(b));
for iter=1:500
    x=D\(b+L*x+U*x);
    error=norm(b-A*x)/norm(b);
    if(error<tol)
        break;
    end
end

function x=nabisect(fname,a,b,e)
%用途：二分法解非线性方程f(x)=0
%格式：x=nabisect(fname,a,b,e).其中，fname为用函数句柄或
% 内嵌函数表达的f(x);a,b为区间端点；e为精度（默认值10-4）；
% x返回解。程序要求函数在两端点值必须异号。
% 中间变量fa,fb,fx引入可以最大限度减少fname调用次数， 从而提高速度。
if nargin<4,e=le-4;end;
fa=fname(a);fb=fname(b);
if fa*fb>0,error('函数在两端点值必须异号');end
x=(a+b)/2
while (b-a)>(2*e)
   fx=fname(x);
   if fa*fx<0,b=x;fb=fx;else a=x;fa=fx;end
   x=(a+b)/2
end
end

function x=nags(A,b,x0,e,N)
%用途:用向量形式(普通储存格式)的高斯-塞德尔迭代解线性方程组Ax=b
%格式:x=nags(A,b,x0,e,N)。A为系数矩阵，b为右端向量，x返回解向量，
%x0为初始向量(默认原点),e为精度(默认le-4),设置迭代次数上限N以防
%发散(默认500)
n=length(b);
if nargin<5,N=500;end
if nargin<4,e=1e-4;end
if nargin<3,x0=zeros(n,1);end
x=x0;x0=x+2*e;
k=0;Al=tril(A);iAl=inv(Al);
while norm(x0-x,inf)>e&&k<N,      
k=k+1;      
x0=x;
x=-iA1*(A-A1)*x0+iA1*b;      
x'
end
if k==N,warning('已达迭代次数上限');end
end

function  [l,u]=nalu(a)
%用途：求可逆方阵的LU分解
%格式：[l,u]=nalu(a)，a为可逆方阵，l为返回的单位下三角矩阵，u为返回的上三角矩阵。
n=length(a);
u=zeros(n,n);
l=eye(n,n);
u(1,:)=a(1,:);
l(2:n,1)=a(2:n,1)/u(1,1);
for k=2:n
u(k,k:n)=a(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);
l(k+1:n,k)=(a(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k))/u(k,k);
end


end

function  a=nalupad(a)
%用途：求可逆方阵的LU分解，紧凑格式
%格式：a=nalupad(a)，a为可逆方阵，l为返回的单位下三角矩阵，u为返回的上三角矩阵。
n=length(a);
a(2:n,1)=a(2:n,1)/u(1,1);
for k=2:n
    a(k,k:n)=a(k,k:n)-a(k,1:k-1)*a(1:k-1,k:n);
    a(k+1:n,k)=(a(k+1:n,k)-a(k+1:n,1:k-1)*a(1:k-1,k))/a(k,k);
end
end

function[g,te]=Newton_f(f,x0,tol)
% x0为迭代初值
% tol为指定误差容限，如果默认为10的-5次方
% 该方法比较实用，希望读者掌握
if(nargin==2)
    tol=1.0e-5
end
df=diff(sym(f));  %计算原函数的导数
x1=x0;
te=0;             %te记次用
w=0.1;            %给定一个误差初值，以便进入循环计算
while(w>tol)
te=te+1;
fx=subs(f,x1);
df=subs(df,x1);
g=x1-fx/df;
w=abs(g-x1);
x1=g;                       %把迭代后的值献给x1
end
end

function xr= secant(fun,x0,x1,D)
 %  secant.m函数为正割法解方程
 %  fun为给定的非线性方程
 %  x0,x1为给定的初始值
 %  D为近似值的误差界
 %  xr为所得的非线性方程的近似解
if nargin ==3;
   D=1e-6;                         %设置默认精度
end
f0 = feval(fun,x0);                %计算x0处的函数值
f1 = feval(fun,x1);                %计算x1处的函数值
while abs(x0-x1)>D;
     x2 = x1- f1*[x0-x1]/[f0-f1];  %正割法迭代公式
     f0 = f1;                      %更新f0数值
     x0 = x1;                      %更新x0数值
     x1 = x2;                      %更新x1数值
     f1 = feval(fun,x1);           %计算x1处的函数值
end
xr = x2;                           %返回根
end

function [x,err,yc]= sfen (f,a,b,eps)
%sfen.m函数用二分法求解非线性方程
%a,b表示求解区间[a,b]的端点，并且满足f(a)*f(b)<0
%f为所求解的非线性方程，可用函数inline生成
%eps表示精度指标；x为所得到的近似解
%err为x的误差估计；yc为函数f在x上的值
ya=feval('f',a);
yb=feval('f',b);
if yb==0
    x=b;
end
if ya*yb>0,
      error('函数在两端点值必须异号');
end
max1=1+round((log(b-a)-log(eps))/log(2));
for k=1:max1
    x=(a+b)/2;
    yc=feval('f',x);
    if yc==0;
        a=x;
        b=x;
        break;
    elseif yb*yc>0
        b=x;
        yb=yc;
    else
        a=x;
        ya=yc;
    end
end
k
x=(a+b)/2
err=abs(b-a)
yc=feval('f',x)
end