计算方法常用计算编程.docx
《计算方法常用计算编程.docx》由会员分享，可在线阅读，更多相关《计算方法常用计算编程.docx（51页珍藏版）》请在冰豆网上搜索。
计算方法常用计算编程
1．向量范数：
clear
x=input（'输入一个向量\nx='）
x1=sum（abs（x））%1-范数
x2=sqrt（sum（abs（x）.^2））%2-范数
xinf=max（abs（x））%无穷-范数
p=input（'多少范数？
输入p='）
xp=（sum（abs（x）.^p））.^（1/p）%p-范数
2．矩阵范数：
clear
a=input（'输入一个矩阵\na='）
x1=max（sum（abs（a）））%1-范数
x2=sqrt（max（eig（a'*a）））%2-范数
xinf=max（sum（abs（a）'））%无穷-范数
xF=sqrt（sum（sum（a.^2）））%F-范数
3．高斯消去法：
设Ax=b,
则
clear
a=input（'输入一个系数矩阵\na='）
b=input（'输入一个自由项\nb='）
n=length（b）;
fork=1:
n-1
fori=k+1:
n
l（i,k）=-a（i,k）/a（k,k）;
forj=k:
n
a（i,j）=a（i,j）+l（i,k）*a（k,j）;
end
b（i）=b（i）+l（i,k）*b（k）;
end
end
zengguang=[a,b]
x（n）=b（n）/a（n,n）;
fork=n-1:
-1:
1
t=0;
forj=k+1:
n
t=t+a（k,j）*x（j）
end
x（k）=（b（k）-t）/a（k,k）;
end
x
输入一个系数矩阵
a=[111;04-1;2-21]
a=
111
04-1
2-21
输入一个自由项
b=[6;5;1]
b=
6
5
1
zengguang=
1116
04-15
00-2-6
x=
123
4．判断一个数字是否是素数。
clear
a=input（'a='）
fori=2:
a-1
ifrem（a,i）==0
sprintf（'%d不是素数',a）
break
end
end
ifa==i+1
sprintf（'%d是素数',a）
end
clear
clc
a=input（'输入系数矩阵a='）
b=input（'输入自由项b='）
n=length（b）;
D=diag（diag（a））
L=zeros（n）;
U=zeros（n）;
fori=1:
n-1
L（i+1:
n,i）=-a（i+1:
n,i）
end
fori=1:
n-1
U（i,i+1:
n）=-a（i,i+1:
n）
end
x1=input（'迭代初始值x1='）
x1=x1';
derta=input（'输入精度要求derta='）
ifderta<=0derta=0.00001
end
%%%%%%%%%%%%%%%%雅克比迭代法矩阵迭代%
t=1;
while（t>derta）
x2=inv（D）*（L+U）*x1+inv（D）*b;
t=max（abs（x1-x2））;
x1=x2;
end
input（'回车后输出雅克比迭代法矩阵迭代的解x='）
x2
%%%%%%%%%%%%%%雅克比迭代法元素迭代法%
t=1;
while（t>derta）
fori=1:
n
tt=0;
forj=1:
i-1
tt=tt+a（i,j）*x1（j）;
end
ttt=0;
forj=i+1:
n
ttt=ttt+a（i,j）*x1（j）;
end
x2（i）=（b（i）-tt-ttt）/a（i,i）;
end
t=max（abs（x1-x2））;
x1=x2;
end
input（'回车后输出雅克比迭代法元素迭代的解x='）
x2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%高斯迭代法矩阵迭代%
t=1;
while（t>derta）
x2=inv（D-L）*U*x1+inv（D-L）*b;
t=max（abs（x1-x2））;
x1=x2;
end
input（'回车后输出高斯迭代法矩阵迭代的解x='）
x2
%%%%%%%%%%%%%%高斯迭代法元素迭代法%
t=1;
while（t>derta）
fori=1:
n
tt=0;
forj=1:
i-1
tt=tt+a（i,j）*x2（j）;
end
ttt=0;
forj=i+1:
n
ttt=ttt+a（i,j）*x1（j）;
end
x2（i）=（b（i）-tt-ttt）/a（i,i）;
end
t=max（abs（x1-x2））;
x1=x2;
end
input（'回车后输出高斯迭代法元素迭代的解x='）
x2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%逐次超松弛迭代法矩阵迭代%
t=1;
w=input（'输入超松弛因子w='）
ifw>0&w<=2
w=w;
elsew=0.15;
end
while（t>derta）
x2=inv（D-w*L）*（（1-w）*D+w*U）*x1+w*inv（D-w*L）*b;
t=max（abs（x1-x2））;
x1=x2;
end
input（'回车后输出超松弛迭代法矩阵迭代的解x='）
x2
%%%%%%%%%%%%%%超松弛迭代法元素迭代法%
t=1;
while（t>derta）
fori=1:
n
tt=0;
forj=1:
i-1
tt=tt+a（i,j）*x2（j）;
end
ttt=0;
forj=i+1:
n
ttt=ttt+a（i,j）*x1（j）;
end
x2（i）=w*（b（i）-tt-ttt）/a（i,i）+x1（i）;
end
t=max（abs（x1-x2））;
x1=x2;
end
input（'回车后输出超松弛迭代法元素迭代的解x='）
x2
%输入系数矩阵a=[-4111;1-411;11-41;111-4]
%输入自由项b=[1;1;1;1]
%迭代初始值x1=[0000]
%输入精度要求derta=0.000046
clear
a=input（'任意输入一个数a='）
fork=2:
a-1
ifrem（a,k）==0
sprintf（'%d不是素数',a）
break
end
end
ifk==a-1
sprintf（'%d是素数',a）
end
%%%向量范数通用程序
clear
x=input（'输入一个向量\nx='）
x1=sum（abs（x））%1-范数
x2=sqrt（sum（abs（x）.^2））%2-范数
xinf=max（abs（x））%无穷-范数
p=input（'多少范数？
输入p='）
xp=（sum（abs（x）.^p））.^（1/p）%p-范数
%%%矩阵范数通用程序
clear
a=input（'输入一个矩阵\na='）
x1=max（sum（abs（a）））%1-范数
x2=sqrt（max（eig（a'*a）））%2-范数
xinf=max（sum（abs（a）'））%无穷-范数
xF=sqrt（sum（sum（a.^2）））%F-范数
%%%高斯顺序消去法通用程序
clear
a=input（'输入一个系数矩阵\na='）
b=input（'输入一个自由项\nb='）
n=length（b）;
fork=1:
n-1
fori=k+1:
n
l（i,k）=-a（i,k）/a（k,k）;
forj=k:
n
a（i,j）=a（i,j）+l（i,k）*a（k,j）;
end
b（i）=b（i）+l（i,k）*b（k）;
end
end
zengguang=[a,b]
x（n）=b（n）/a（n,n）;
fork=n-1:
-1:
1
t=0;
forj=k+1:
n
t=t+a（k,j）*x（j）
end
x（k）=（b（k）-t）/a（k,k）;
end
x
clear
clc
x0=input（'输入初值x0='）
ifabs（x0）>0
x0=x0;
elsex0=1.5
end
x=sym（'x'）;
%f='x-x^3-4*x^2+10'
f=input（'输入迭代函数f='）
t=1;
pp=1;
%%%%%%下面判断迭代是否收敛
g=diff（f）;
x=x0;
gg=eval（g）
ifabs（gg）>1
x1='迭代格式不收敛'
else
derta=input（'输入迭代精度derta='）
ifabs（derta）>0
derta=derta;
elsederta=0.000001;
end
t=1;
x=x0;
while（t>derta）
pp=pp+1;
x1=eval（f）;
t=abs（x1-x）;
x=x1;
ifpp>10000%如果不收敛，则迭代10000次后自动终止
x1='迭代超过10000次';
break;
end
end
end
p
x1
%f='x-x^3-4*x^2+10'
%f='1/2*（10-x^3）^（1/2）'
%f='（10/x-4*x）^（1/2）'
%f='（10/（4+x））^（1/2）'
%f='x-（x^3+4*x^2-10）/（3*x^2+8*x）'
%%%%%%%%%%%%%%%%结果展示
%输入初值x0=1.5
%x0=
%1.5000
%输入迭代函数f='x-x^3-4*x^2+10'
%f=
%x-x^3-4*x^2+10
%输入迭代精度derta=0.00002
%derta=
%2.0000e-005
%x1=
%NaN
clear
clc
x0=input（'输入初值x0='）
ifabs（x0）>0
x0=x0;
elsex0=1.5
end
x=sym（'x'）;
%f='x-x^3-4*x^2+10'
f=input（'输入迭代函数f='）
t=1;
pp=1;
%%%%%%下面判断迭代是否收敛
g=diff（f）;
x=x0;
gg=eval（g）;
ifabs（gg）>1
x1='迭代格式不收敛'
else
derta=input（'输入迭代精度derta='）
ifabs（derta）>0
derta=derta;
elsederta=0.000001;
end
t=1;
x=x0;
w=1/（1-gg）;%加权因子
while（t>derta）
pp=pp+1;
x1=w*x+（1-w）*eval（f）;%两个组合的叠加
t=abs（x1-x）;
x=x1;
ifpp>10000%如果不收敛，则迭代10000次后自动终止
x1='迭代超过10000次';
break;
end
end
end
p
x1
%f='x-x^3-4*x^2+10'
%f='1/2*（10-x^3）^（1/2）'
%f='（10/x-4*x）^（1/2）'
%f='（10/（4+x））^（1/2）'
%f='x-（x^3+4*x^2-10）/（3*x^2+8*x）'
%%%%%%%%%%%%%%%%
clear
clc
x0=input（'输入初值x0='）
ifabs（x0）>0
x0=x0;
elsex0=1.5
end
x=sym（'x'）;
%f='x-x^3-4*x^2+10'
f=input（'输入迭代函数f='）
t=1;
pp=1;
%%%%%%下面判断迭代是否收敛
g=diff（f）;
x=x0;
gg=eval（g）;
ifabs（gg）>1
x1='迭代格式不收敛'
else
derta=input（'输入迭代精度derta='）
ifabs（derta）>0
derta=derta;
elsederta=0.000001;
end
t=1;
x=x0;
w=1/（1-gg）;%加权因子
while（t>derta）
pp=pp+1;
x2=x;
y=eval（f）;
x=y;
z=eval（f）;
x=x2;
x1=x-（y-x）^2/（z-2*y+x）;%两个组合的叠加
t=abs（x1-x）;
x=x1;
ifpp>10000%如果不收敛，则迭代10000次后自动终止
x1='迭代超过10000次';
break;
end
end
end
p
x1
%f='x-x^3-4*x^2+10'
%f='1/2*（10-x^3）^（1/2）'
%f='（10/x-4*x）^（1/2）'
%f='（10/（4+x））^（1/2）'
%f='x-（x^3+4*x^2-10）/（3*x^2+8*x）'
%%%%%%%%%%%%%%%%
%乘幂法
clear
clc
A=input（'输入一个n阶方阵A='）
n=length（A（1,:
））
v0=ones（1,n）';
u0=v0/max（v0）
q=1;
while（q>0.00000001）
v1=A*u0
j=1;
t=1;
fori=1:
n
ifabs（v1（i））>abs（v1（j））j=i;
end
end
fori=1:
n
ifabs（v0（i））>abs（v0（t））t=i;
end
end
u1=v1./v1（j）
q=abs（abs（v1（j））-abs（v0（t）））;
u0=u1;
v0=v1;
end
j=1;
fori=1:
n
ifabs（v1（i））>abs（v1（j））j=i;
end
end
lamda=v1（j）
tezhengxiangliang=v1'./v1（j）
%输入一个n阶方阵A=[-1233;31-2;3-27]
clear
clc
x0=input（'输入起始节点坐标x0='）
h=input（'输入步长h='）
y=input（'输入节点坐标函数值f（x）='）
x2=input（'输入所要计算的节点x2='）
symst
n=length（y）;
fori=1:
n
x1（i）=x0+h*（i-1）;
end
%%%%%%%%%%%%%%%%%%%%差分的求法
fori=2:
n
f1（i,1）=（y（i）-y（i-1））;
end
fori=2:
n
forj=i+1:
n
f1（j,i）=（f1（j,i-1）-f1（j-1,i-1））;
end
end
f1=[y',f1]%输出带0阶差分的差分表格
%%%%%%%%%%%%%%%%%%%%%%%%%%等距节点插值函数中牛顿向前插值函数
Newton1=f1（1,1）;
ttt=1;
fori=2:
n
tt=1;
forj=1:
i-1
tt=tt*（t-j+1）;
end
ttt=ttt*（i-1）;
Newton1=Newton1+f1（i,i）*tt/（ttt）
end
fprintf（'等距节点插值函数中牛顿向前插值函数为\n'）
expand（Newton1）
t=（x2-x0）/h;
p=eval（Newton1）;
fprintf（'等距节点插值函数中牛顿向前插值函数为在所求点x2的函数值为\n'）
p
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%输入起始节点坐标x0=0.4
%输入步长h=0.1
%输入节点坐标函数值f（x）=[0.389420.479430.56464]
%输入所要计算的节点x2=0.43251
%Newton插值函数为为:
1+17/6*x-2/3*x^2+1/6*x^4-1/3*x^3
%Newton插值函数在所求点x2的函数值为:
23
clear
clc
x0=input（'输入起始节点坐标x0='）
h=input（'输入步长h='）
y=input（'输入节点坐标函数值f（x）='）
x2=input（'输入所要计算的节点x2='）
symst
n=length（y）;
fori=1:
n
x1（i）=x0+h*（i-1）;
end
%%%%%%%%%%%%%%%%%%%%差分的求法
fori=2:
n
f1（i,1）=（y（i）-y（i-1））;
end
fori=2:
n
forj=i+1:
n
f1（j,i）=（f1（j,i-1）-f1（j-1,i-1））;
end
end
f1=[y',f1]%输出带0阶差分的差分表格
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%等距节点插值函数中牛顿向后插值函数
Newton2=f1（n,1）;
ttt=1;
fori=2:
n
tt=1;
forj=1:
i-1
tt=tt*（t-j+1）;
end
ttt=ttt*（i-1）;
Newton2=Newton2+f1（n,i）*tt/（ttt）
end
fprintf（'等距节点插值函数中牛顿向后插值函数为\n'）
expand（Newton2）
t=（x2-x0）/h;
p=eval（Newton2）;
fprintf（'等距节点插值函数中牛顿向后插值函数为在所求点x2的函数值为\n'）
p
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%输入起始节点坐标x0=0.4
%输入步长h=0.1
%输入节点坐标函数值f（x）=[0.389420.479430.56464]
%输入所要计算的节点x2=0.43251
%Newton插值函数为为:
1+17/6*x-2/3*x^2+1/6*x^4-1/3*x^3
%Newton插值函数在所求点x2的函数值为:
23
clear
clc
formatlong
symsx
f=input（'输入被积函数f='）
a=input（'输入积分区间[ab]='）
%%%%%%%%%%%%%%%%%%%%%梯形积分公式
（'梯形积分公式'）
x=a
（1）;
f1=eval（f）;
x=a
（2）;
f2=eval（f）;
fun=（f1+f2）*（a
（2）-a
（1））/2
%%%%%%%%%%%%%%%%%%%%%辛普森积分公式
（'辛普森积分公式'）
x=a
（1）;
f1=eval（f）;
x=（a
（1）+a
（2））/2;
f2=eval（f）;
x=a
（2）;
f3=eval（f）;
fun=（f1+4*f2+f3）*（a
（2）-a
（1））/6
%%%%%%%%%%%%%%%%%%%%%科特斯积分公式
（'科特斯积分公式'）
h=（a
（2）-a
（1））/4;
x=a
（1）;
f1=eval（f）;
x=a
（1）+h;
f2=eval（f）;
x=a
（1）+2*h;
f3=eval（f）;
x=a
（1）+3*h;
f4=eval（f）;
x=a
（1）+4*h;
f5=eval（f）;
fun=（7*f1+32*f2+12*f3+32*f4+7*f5）*（a
（2）-a
（1））/90
%实例验证p176例7.5
%输入被积函数f=sqrt（x）
%输入积分区间[ab]=[0.51]
clear
clc
a=input（'输入系数矩阵a='）
b=input（'输入自由项b='）
x1=input（'输入初始值x1='）
n=length（b）;
ifsum（abs（x1））>0
x1=x1;
elsex1=zeros（1,n）;
end
t=2;
%%%%%%%%%%%%%%%%%雅克比迭代法
while（t>0.00001）
fori=1:
n
q=0;
forj=1:
i-1
q=q+a（i,j）*x1（j）;
end
p=0;
forj=i+1:
n
p=p+a（i,j）*x1（j）;
end
x2（i）=（b（i）-q-p）/a（i,i）;
end
t=max（abs（x2-x1））;
x1=x2;
end
x2
%%%%%%%%%%%%%%%%%%%%%%%%高斯赛德尔迭代法
while（t>0.00001）
fori=1:
n
q=0;
forj=1:
i-1
q=q+a（i,j）*x2（j）;
end
p=0;
forj=i+1:
n
p=p+a（i,j）*x1（j）;
end
x2（i）=（b（i）-q-p）/a（i,i）;
end
t=max（abs（x2-x1））;
x1=x2;
end
x2
%%%%%%%%%%%%%%%%%%%%%%%%%%%超松弛迭代法
w=input（'输入松弛因子（0ifw>0w=w;
elsew=0.5
end
while（t>0.00001）
fori=1:
n
q=0;
forj=1:
i-1
q=q+a（i,j）*x2（j）;
end
p=0;
forj=i+1:
n
p=p+a（i,j）*x1（j）;
end
x2（i）=w*（b（i）-q-p）/a（i,i）+x1（i）;
end
t=max（abs（x2-x1））;
x1=x2;
end
x2
%%%%%%%%%%%
%输入系数矩阵a=[-4111;1-411;11-41;111-4]
%a=
%-4111
%1-411
%11-41
%111-4
%输入自由项b=[1111]
%b=
%1111
%输入初始值x1=[0000]
%x1=
%0000
%雅克比迭代结果
%x2=
%-1.0000-1.0000-1.0000-1.0000
%高斯赛德尔迭代结果
%x2=
%-1.0000-1.0000-1.0000-1.0000
%
升级会员 