北科大数值分析考试程序.docx
《北科大数值分析考试程序.docx》由会员分享，可在线阅读，更多相关《北科大数值分析考试程序.docx（28页珍藏版）》请在冰豆网上搜索。
北科大数值分析考试程序
%牛顿法
function[gen,time]=Newton（f,x0,tol）
if（nargin==2）
tol=1.0e-5;
end
%计算原函数的导数
df=diff（sym（f））;
x1=x0;
%time记次用
time=0;
wucha=0.1;%给定一个误差初值以方进入循环计算
while（wucha>tol）
time=time+1;
fx=subs（f,x1）;
df=subs（df,x1）;
gen=x1-fx/df;
wucha=abs（gen-x1）;
x1=gen;%把迭代后的值赋给x1
end
end
%重根的加速迭代问题1
function[gen,time]=Newton_Tay（f,x0,tol）
%x0迭代初值；tol为指定误差容限，若缺省默认为10的-5次方
if（nargin==2）
tol=1.0e-5;
end
%计算原函数的导数
df=diff（sym（f））;
df2=diff（df）;%2阶导数
x1=x0;
%time记次数
time=0;
error=0.1;%给定一个误差初值，以方便进入循环计算
while（error>tol）
time=time+1;
fx=subs（f,x1）;
df=subs（df,x1）;
df2=subs（df2,x1）;
gen=x1-fx*df/（df^2-fx*df2）;%迭代公式
error=abs（gen-x1）;
x1=gen;%把迭代后的值赋给x1
end
end
%割线法计算非线性方程
function[gen,time]=Secant（f,a,b,tol）
if（nargin==3）
tol=1.0e-5;
end
time=0;
error=0.1;
fa=subs（sym（f）,a）;
fb=subs（sym（f）,b）;
gen=a-（b-a）*fa/（fb-fa）;
while（error>tol）
time=time+1;
x1=gen;
fx=subs（sym（f）,x1）;
s=fx*fa;
if（s>0）
gen=b-（x1-b）*fb/（fx-fb）;
else
gen=a-（x1-a）*fa/（fx-fa）;
end
error=abs（gen-x1）;
end
end
%Steffensen’smethod加速其收敛
function[gen,time]=Steff（fun,x0,tol）
if（nargin==2）
tol=1.0e-5;
end
%设置误差初值
time=0;
error=0.1;
gen=x0;
while（error>tol）
x1=gen;
y=subs（fun,x1）+x1;
z=subs（fun,y）+y;
%加速公式
if（z-2*y+x1==0）
gen=z;
else
gen=x1-（y-x1）^2/（z-2*y+x1）;
end
error=abs（gen-x1）;
time=time+1;%迭代加一次的记录
end
gen;%计算结果
time;%迭代次数
%求重根g（x）=x-m*f/f'
%Newtonmethod
function[gen,time]=chonggen2（f,x0,m,tol）
if（nargin==2）
tol=1.0e-5;
end
%计算原函数的导数
df=diff（sym（f））;
x1=x0;
time=0;
error=0.1;
while（error>tol）
time=time+1;
fx=subs（f,x1）;
df=subs（df,x1）;
gen=x1-m*fx/df;
error=abs（gen-x1）;
x1=gen;
end
end
%二分法
functiongen=Bisection（f,a,b,tol）
%f为方程f（x）=0中的f（x）
%如果输入变量缺省则默认误差为1E-3
if（nargin==3）
tol=1.0e-3;
end
gen=compute_gen（f,a,b,tol）;
functionr=compute_gen（f,a,b,tol）
%计算左端点函数值
fa=subs（f,a）;
%右端函数值
%fb=subs（f,b）;
%区间中点函数值
fzd=subs（f,（a+b）/2）;
%sub函数R=subs（S,old,new）
%其中S为符号表达式，old为老变量，new为新变量
%sub把老变量替换为新变量
if（fa*fzd>0）
t=（a+b）/2;
%采用递归方法
r=compute_gen（f,t,b,tol）;
else
if（fa*fzd==0）
r=（a+b）/2;
else
if（abs（b-a）<=tol）
r=（b+3*a）/4;
else
s=（a+b）/2;
%结果
r=compute_gen（f,a,s,tol）;
end
end
end
%不动点迭代（Picard迭代）
function[x,time]=Picard（f,x0,tol）
%结果给出迭代次数
%x0为迭代初值
%tol为误差容限
if（nargin==2）
tol=1.0e-5;
end
%缺省情况下误差容限为十的负五次方
wucha=0.5;%设置误差初值
x1=x0;%x1与x0为前后两次计算结果
time=0;%用于记录迭代次数
while（wucha>tol）
x1=subs（f,x0）;
%迭代计算
wucha=abs（x1-x0）;
x0=x1;%更新x0的值在循环中这一句非常重要
time=time+1;
%记下迭代次数
end
x=x1;
%methodofFalsePosition
function[root,time]=Position（f,a,b,tol）
formatlong;%formatlong显示15位双精度
if（nargin==3）
tol=1.0e-4;
end
time=0;
f1=subs（sym（f）,findsym（sym（f））,a）;
f2=subs（sym（f）,findsym（sym（f））,b）;
%sym则是将字符或者数字转换为字符
%findsym（s,n）可以找出表达式s中n个与x接近的变量
%R=subs（S,old,new）利用new的值代替符号表达式中old的值，old为符号变量或是字符串变量名。
if（f1==0）
root=a;
end
if（f2==0）
root=b;
end
if（f1*f2>0）
disp（'两端点的函数值乘积大于0！
'）;
return;
else
tol1=1;
r1=a;%迭代初始值
r2=b;
fv=subs（sym（f）,findsym（sym（f））,a）;
while（tol1>tol）
time=time+1;
f2=subs（sym（f）,findsym（sym（f））,r2）;
root=r2-（r2-r1）*f2/（f2-fv）;
fr=subs（sym（f）,findsym（sym（f））,root）;
if（f2*fr<0）
tol1=abs（root-r2）;
r1=r2;
r2=root;
fv=subs（sym（f）,findsym（sym（f））,r1）;
else
tol1=abs（root-r2）;
r2=root;
fv=0.5*subs（sym（f）,findsym（sym（f））,r1）;
end
end
end
formatshort;
%x0为插值多项式在此点的近似，可取不同的值
function[Q,w]=Neville（x0,n,x,f）
ifn==5
fx=[subs（f,x
（1））,subs（f,x
（2））,subs（f,x（3））,subs（f,x（4））,subs（f,x（5））,subs（f,x（6））]
fori=1:
6
Q（i,1）=fx（i）;
end
end
ifn==10
fx=[subs（f,x
（1））,subs（f,x
（2））,subs（f,x（3））,subs（f,x（4））,subs（f,x（5））,subs（f,x（6））,subs（f,x（7））,subs（f,x（8））,subs（f,x（9））,subs（f,x（10））,subs（f,x（11））]
fori=1:
11
Q（i,1）=fx（i）;
end
end
k=1;
forx0=-1:
0.01:
1
fori=1:
n
forj=1:
i
Q（i+1,j+1）=（（x0-x（i-j+1））*Q（i+1,j）-（x0-x（i+1））*Q（i,j））/（x（i+1）-x（i-j+1））;
end
end
w（k）=Q（i+1,j+1）;
k=k+1;
end
%lagrange插值方法
functions=Lagrange（x,y,t）
%采用符号推导，这样可以给出插值具体公式
symsp;
%读取x向量维数
n=length（x）;
s=0;
for（k=1:
n）
la=y（k）;
%构造基函数
for（j=1:
k-1）
la=la*（p-x（j））/（x（k）-x（j））;
end;
for（j=k+1:
n）
la=la*（p-x（j））/（x（k）-x（j））;
end;
s=s+la;
simplify（s）;%simplify（S）:
对表达式S进行化简。
end
%对输入参数的个数做判断，若只有两个参数，直接给出插值多项式
%如果三个参数则给出插值点的插值结果
%第三个参数可以为向量
if（nargin==2）
s=subs（s,'p','x'）;
%展开多项式
s=collect（s）;
%把系数取到6位精度表达
s=vpa（s,4）;%变量精度计算
else
%读取t长度
m=length（t）;
%分别对t的每一个分量插值
fori=1:
m
temp（i）=subs（s,'p',t（i））;
end
%得到的是系列插值点的插值结果
%既得到的是向量，赋值给s
s=temp;
end
%牛顿插值方法
%如果参数要计算的点列参数缺省，则直接给出插值多项式
%其中要计算的插值点可以是向量
%就是说可以对一系列点进行一次计算完成
functions=newtoninter（x,y,t）
%定义符号变量
symsp;
s=y
（1）;
coefficient=0;
dxs=1;
%读取x的长度
n=length（x）;
%构造Newtoninterpolationmethod
for（i=1:
n-1）
for（j=i+1:
n）
coefficient（j）=（y（j）-y（i））/（x（j）-x（i））;
end
temp1（i）=coefficient（i+1）;
dxs=dxs*（p-x（i））;
s=s+temp1（i）*dxs;
y=coefficient;
end
simplify（s）;
%如果插值点函数缺省则直接输出多项式
if（nargin==2）
s=subs（s,'p','x'）;
s=collect（s）;
s=vpa（s,4）;
else
%读取要插值点向量长度
%可以直接对多点插值计算
m=length（t）;
fori=1:
m
temp（i）=subs（s,'p',t（i））;
end
%得到的是系列插值点的插值结果
%既得到的是向量，赋值给s
s=temp;
end
%Romberg数值求解积分
function[R,k,T]=Romb（fun,a,b,tol）
%fun被积函数
%a,b积分限
%R积分值
%k迭代次数，T整个迭代过程
k=0;
n=1;
h=b-a;
T=h/2*（Romb_f（a）+Romb_f（b））;
error=1;
whileerror>=tol
k=k+1;
h=h/2;
tmp=0;
fori=1:
n
tmp=tmp+Romb_f（a+（2*i-1）*h）;
end
T（k+1,1）=T（k）/2+h*tmp;
forj=1:
k
T（k+1,j+1）=T（k+1,j）+（T（k+1,j）-T（k,j））/（4^j-1）;
end
n=n*2;
error=abs（T（k+1,k+1）-T（k,k））;
end
R=T（k+1,k+1）;
function[f,k,t]=Romb_f（x）
f=sqrt（1+（cos（x））.^2）;
[R,k,T]=Romb（'sqrt（1+（cos（x））.^2）',0,48,1.e-8）
x=0:
0.02:
50;
y=sqrt（1+（cos（x））.^2）;
area（x,y）
grid
title（'Romberg积分'）
%Romberg积分方法
%a,b为积分区间，tol为误差容限
%s为最后的积分面积
function[s,k]=Romberg（a,b,tol）
n=1;
h=b-a;
%设置设计误差初值
delt=1;
x=a;
k=0;
R=zeros（4,4）;
R（1,1）=h*（Rombg_f（a）+Rombg_f（b））/2;
whiledelt>tol
%如果两次计算的差值大于给定误差则进入循环
k=k+1;
h=h/2;
s=0;
forj=1:
n
x=a+h*（2*j-1）;
s=s+Rombg_f（x）;
end
R（k+1,1）=R（k,1）/2+h*s;
n=2*n;
fori=1:
k
R（k+1,i+1）=（（4^i）*R（k+1,i）-R（k,i））/（4^i-1）;
end
%前后两次值的差
delt=abs（R（k+1,k）-R（k+1,k+1））;
end
s=R（k+1,k+1）;
k=k+1;
%Romberg方法试验函数
function[f,t]=Rombg_f（x）
f=sqrt（1+（cos（x））.^2）;
%Romberg积分的主函数
[s,t]=Romberg（0,48,1.0e-6）
x=0:
0.02:
48;
y=sqrt（1+（cos（x））.^2）;
area（x,y）
grid
%Eulermethod
function[t,x]=Euler（fun,t0,t1,x0,h）
%fun表示f（t,x）
%t0,t1表示自变量的初值和终值
%x0表示函数在x0处的值
N=（t1-t0）/h;
t=t0+[0:
N]'*h;%t变化
x（1,:
）=x0';%赋初值
fork=1:
N
f=feval（fun,t（k）,x（k,:
））;%将t，x带入fun函数中
f=f';
x（k+1,:
）=x（k,:
）+h*f;
end
functionf=Euler_fun（t,x）
f=-x+t+1;
%ModifiedEulermethod
function[t,x]=ModifiedEuler（fun,t0,tt,x0,h）
N=（tt-t0）/h;
t=t0+[0:
N]'*h;
x（1,:
）=x0';
fori=1:
N
f1=h*feval（fun,t（i）,x（i,:
））;
f1=f1';
f2=h*feval（fun,t（i+1）,x（i,:
）+f1）;
f2=f2';
x（i+1,:
）=x（i,:
）+1/2*（f1+f2）;
end
functionf=ModifiedEuler（t,x）
f=-x+t+1;
%Runge-Kutta（orderfour）methodforsystemsofDifferentialequations
functionRK（A,x,h,y0）
%A为函数组
%x为轴上的各点
%y0为初值
i=1;
y（i,:
）=y0;
m=length（A）;
b=x（length（x））;
whilex（i）a=[x（i）,y（i,:
）];
forl=1:
m
k1（l）=eval（A{l},a）;%eval（）函数的功能就是将括号内的字符串视为语句并运行
end
a=[x（i）+h/2,y（i,:
）+h/2*k1];
forl=1:
m
k2（l）=eval（A{l},a）;
end
a=[x（i）+h/2,y（i,:
）+h/2*k2];
forl=1:
m
k3（l）=eval（A{l},a）;
end
a=[x（i）+h,y（i,:
）+h*k3];
forl=1:
m
k4（l）=eval（A{l},a）;
end
y（i+1,:
）=y（i,:
）+h/6*（k1+2*k2+2*k3+k4）;
i=i+1;
end
y
%Runge-KuttaMethodtosolvetheODE
function[t,y]=RK4（fun,t0,tt,y0,h,varargin）
%nargin对应varargin的个数,varargin用来存入输入变量
%y0初始时函数值
ifnargin<=4
h=0.001
end%如果没有输入h的话，那么默认计算步N（tt-t0）/0.001
y（1,:
）=y0（:
）';
N=（tt-t0）/h;
t=t0+[0:
N]'*h;
fork=1:
N
%RK4中k1
f1=h*feval（fun,t（k）,y（k,:
））;%将t，y的值一次带入fun函数中
f1=f1（:
）';
%RK4中k2
f2=h*feval（fun,t（k）+h/2,y（k,:
）+f1/2）;
f2=f2（:
）';
%RK4中k3
f3=h*feval（fun,t（k）+h/2,y（k,:
）+f2/2）;
f3=f3（:
）';
%RK4中k4
f4=h*feval（fun,t（k）+h,y（k,:
）+f3）;
f4=f4（:
）';
%结果
y（k+1,:
）=y（k,:
）+（f1+2*f2+2*f3+f4）/6;
end
functionf=RK4_fun（t,y）
f=-y+t+1;
%GaussianeliminationwithPartialPivoting
function[RA,RB,n,X]=Column（A,b）
B=[Ab];
n=length（b）;
RA=rank（A）;
RB=rank（B）;
ifRB-RA>0
disp（'因为RA~=RB，所以此方程组无解。
'）
return
end
ifRA==RB
ifRA==n
disp（'由于RA=RB=n，所以此方程组有唯一解'）
X=zeros（n,1）;
C=zeros（1,n+1）;
forp=1:
n-1
[Y,j]=max（abs（B（p:
n,p）））;
C=B（p,:
）;
B（p,:
）=B（j+p-1,:
）;
B（j+p-1,:
）=C;
fork=p+1:
n
m=B（k,p）/B（p,p）;
B（k,p:
n+1）=B（k,p:
n+1）-m*B（p,p:
n+1）;
end
end
b=B（1:
n,n+1）;
A=B（1:
n,1:
n）;
X（n）=b（n）/A（n,n）;
fort=n-1:
-1:
1
X（t）=（b（t）-sum（A（t,t+1:
n）*X（t+1:
n）））/A（t,t）;
end
else
disp（'由于RA=RB'）
end
end
%GaussianeliminationwithPartialPivoting----按列选主元
A=[1.192.11-1001
14.2-0.12212.2-1
0100-99.91
15.30.110-13.1-1];
b=[1.123.442.154.16]';
[RA,RB,n,X]=Column（A,b）
%Gaussianeliminitation
%A是系数矩阵，B是增广矩阵
function[RA,RB,n,X]=Gauss（A,b）
B=[Ab];
n=length（b）;
RA=rank（A）;
RB=rank（B）;
ifRB-RA>0
disp（'因为RA~=RB，所以此方程组无解。
'）
return
end
ifRA==RB
ifRA==n
disp（'由于RA=RB=n，所以此方程组有唯一解'）
X=zeros（n,1）;
C=zeros（1,n+1）;
forp=1:
n-1
fork=p+1:
n
m=B（k,p）/B（p,p）;
B（k,p:
n+1）=B（k,p:
n+1）-m*B（p,p:
n+1）;
end
end
b=B（1:
n,n+1）;
A=B（1:
n,1:
n）;
X（n）=b（n）/A（n,n）;
fort=n-1:
-1:
1
X（t）=（b（t）-sum（A（t,t+1:
n）*X（t+1:
n）））/A（t,t）;
end
else
disp（'由于RA=RB'）
end
end
%Gausseliminiation解线性方程组
A=[1.192.11-1001
14.2-0.12212.2-1
0100-99.91
15.30.110-13.1-1];
b=[1.123.442.154.16]';
[RA,RB,n,X]=Gauss（A,b）
%xx=[0.176825300.01269269-0.02065405-1.18260870]';
%usingGauss-Seidelmethodtosolvethesystems
function[X,k]=GS（A,b,x0,tol）
max1=300;
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
T=（D-L）\U;
C=（D-L）\b;
X=T*x0+C;
k=1;
whilenorm（X-x0）>=tol
x0=X;
X=T*x0+C
升级会员 