版权声明:转载请注明出处:https://blog.csdn.net/qq1013459920 https://blog.csdn.net/qq1013459920/article/details/84643272
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