计算方法matlab上机报告流程图源代码等.docx
《计算方法matlab上机报告流程图源代码等.docx》由会员分享，可在线阅读，更多相关《计算方法matlab上机报告流程图源代码等.docx（12页珍藏版）》请在冰豆网上搜索。
计算方法matlab上机报告流程图源代码等
上机报告
1.共轭梯度法
（1）计算过程如下：
第一步：
取初始向量
计算
第
步：
计算
（2）
symsn;%定义一个变量n
A=input（'请输入矩阵A'）
b=input（'请输入矩阵b'）
n=size（A）
X=zeros（n）;
D=zeros（n）;
R=zeros（n）;
x=X（:
1）;%将矩阵X的第一列赋值给x作为初始向量
r=b-A*x;%将x代入求得初始非零残向量
R（:
1）=r;
d=r;%求得初始方向向量
D（:
1）=d;
fori=1:
n;%利用循环进行迭代求得各向量
AF（i）=R（:
i）'*R（:
i）/（D（:
i）'*A*D（:
i））;
X（:
i+1）=X（:
i）+AF（i）*D（:
i）;
R（:
i+1）=b-A*X（:
i+1）;
BT（i）=norm（R（:
i+1））^2/norm（R（:
i））^2;
D（:
i+1）=R（:
i+1）+BT（i）*D（:
i）;
end
x=X（:
i+1）%输出计算结果
2.龙贝格积分
s=input（'请输入被积函数表达式：
f（x）=','s'）%键盘输入被积函数表达式
f=inline（s）
a=input（'请输入积分下限a='）%输入积分下限a
b=input（'请输入积分上限b='）%输入积分上限b
f0=f（a）;%求下限值f（a）
f1=f（b）;%求上限值f（b）
j=zeros（8,4）;%定义一个矩阵用来存放T,S,C,R值
j（1,1）=（b-a）*0.5*（f0+f1）;%计算出T1
fori=2:
8;%对i循环赋值
t=2^（i-2）;
fj=zeros（t,1）;%构造一个矩阵为接下来存放f（a+（2i-1）*（b-a）/2^（k+1））的值做准备
form=1:
t;
fj（m,1）=f（a+（2*m-1）*（b-a）/2^（i-1））;%将计算得到的f（a+（2i-1）*（b-a）/2^（k+1））的值赋值给对应的矩阵元素
ff=sum（fj）;%对矩阵所有元素求和
end
j（i,1）=0.5*j（i-1,1）+（1/2^（i-1））*ff;%得到所有的T值
j（i,2）=j（i,1）+（1/3）*（j（i,1）-j（i-1,1））;%得到所有的S值
j（i,3）=j（i,2）+（1/3）*（j（i,2）-j（i-1,2））;%得到所有的C值
j（i,4）=j（i,3）+（1/3）*（j（i,3）-j（i-1,3））;%得到所有的R值
s=j（i,4）-j（i-1,4）;
end
j%输出计算结果表
If=vpa（j（i,4）,7）%得到精确的积分结果
3.三样条插值
X=zeros（1,6）;%定义一些下面将用到的矩阵
y=zeros（1,6）;
f=zeros（6）;
l=zeros（1,6）;
fori=1:
1:
6
X（i）=（（i-1）*2）/5-1;%将区间等分取点
y（i）=1/（1+25*X（i）*X（i））;%得到对应点的函数值
end
forj=2:
6
f（1,1）=y
（1）;
f（j,1）=y（j）;
fork=2:
j
f（j,k）=（f（j,k-1）-f（j-1,k-1））/（X（j）-X（j-k+1））;%利用循环求差商
end
end
f%得到差商表，其中对角线上的为需要的差商值
symsx;
l=[x,x,x,x,x,x];
g=l-X;
n=g;
fort=2:
6
n（1,t）=n（1,t-1）*n（1,t）;
end
N=zeros（1,1）;
N=f（1,1）;
forr=1:
5
N=N+n（1,r）*f（r+1,r+1）;%利用循环求5次牛顿插值多项式
end
N%得到5次牛顿插值多项式
P=zeros（1,101）;
X1=[-1:
0.02:
1];%取间距为0.02的等分点作为作图的横坐标
foru=1:
101
x=X1（u）;
subs（N）;
p（1,u）=subs（N（1,1））;%用上面所求出的5次牛顿插值多项式算得所取等分点的函数值作为纵坐标
end
plot（X1,p,'-c'）%绘制5次牛顿插值曲线
holdon
XX=zeros（1,11）;%定义一些下面将用到的矩阵
yy=zeros（1,11）;
F=zeros（11）;
L=zeros（1,11）;
forI=1:
1:
11
XX（I）=（（I-1）*2）/10-1;%将区间等分取点
yy（I）=1/（1+25*XX（I）*XX（I））;%得到对应点的函数值
end
forJ=2:
11
F（1,1）=yy
（1）;
F（J,1）=yy（J）;
forK=2:
J
F（J,K）=（F（J,K-1）-F（J-1,K-1））/（XX（J）-XX（J-K+1））;%利用循环求差商
end
end
F%得到差商表，其中对角线上的为需要用的差商值
symsx1;
L=[x1,x1,x1,x1,x1,x1,x1,x1,x1,x1,x1];
G=L-XX;
M=G;
forT=2:
11
M（1,T）=M（1,T-1）*M（1,T）;
end
H=F（1,1）;
forR=1:
10
H=H+M（1,R）*F（R+1,R+1）;%利用循环求10次牛顿插值多项式
end
H%得到10次牛顿插值多项式
P=zeros（1,101）;
XX1=[-1:
0.02:
1];%取间距为0.02的等分点作为作图的横坐标
foru1=1:
101
x1=XX1（u1）;
subs（H）;
P（1,u1）=subs（H）;%用上面所求出的10次牛顿插值多项式算得所取等分点的函数值作为纵坐标
end
plot（XX1,P,'-g'）
holdon%绘制10次牛顿插值曲线
x3=[-1:
0.02:
1];%取间距为0.02的等分点作为作图的横坐标
y3=1./（1+25*x3.^2）;%用原函数计算得到所取等分点的纵坐标
plot（x3,y3,'-k'）
holdon%绘制原函数f（x）的曲线
h=0.2;
u=0.5;
z=0.5;
d=zeros（11,1）;
d
（1）=6*（（yy
（2）-yy
（1））/0.2-50/（26*26））/0.2;
d（11）=6*（-50/（26*26）-（yy（11）-yy（10））/0.2）/0.2;
a=0.5*ones（1,10）;
b=2*ones（1,11）;
c=diag（a,1）+diag（a,-1）+diag（b）;
c（1,2）=1;
c（11,10）=1;
e=zeros（11,1）;%得到三次样条插值第二种边界条件对应的严格三对角占优矩阵
e（1,1）=d
（1）;
e（11,1）=d（11）;
forii=2:
10
e（ii,1）=F（ii+1,3）;
end
e;%得到d0,d1,….dn用矩阵e表示
ff=horzcat（c,e）%利用d0,d1,…..dn与上面所得的三对角占优矩阵构造增广矩阵ff
y0=zeros（11,1）;%定义一些下面将要用的矩阵
uu=zeros（11,12）;
uu（1,1）=ff（1,1）;
ll=zeros（11,11）;
ll（1,1）=1;
y0（1,1）=ff（1,12）;
fortt=2:
11
ll（tt,tt）=1;
ll（tt,1）=ff（tt,1）/uu（1,1）;
l（tt,tt-1）=ff（tt,tt-1）/uu（tt-1,tt-1）;
forjj=tt:
12
uu（1,jj）=ff（1,jj）;
uu（tt,jj）=ff（tt,jj）-ll（tt,tt-1）*uu（tt-1,jj）;%对增广矩阵进行LU分解，同时进行追法求y0
y0（tt,1）=uu（tt,12）;
end
end
xx=zeros（11,1）;
forkk=10:
-1:
1
xx（11,1）=y0（11,1）/uu（11,12）;
xx（kk,1）=（y0（kk,1）-uu（kk,kk+1）*xx（kk+1,1））/uu（kk,kk）;%用赶法求得xx
end
M0=xx;
Sx=zeros（10,1）;
fori0=1:
10
symsx2;
Sx=sym（Sx）;
Sx（i0）=（（XX（i0+1）-x2）*（XX（i0+1）-x2）*（XX（i0+1）-x2）/（6*h））*M0（i0）+...
（（x2-XX（i0））*（x2-XX（i0））*（x2-XX（i0））/（6*h））*M0（i0+1）+...
（yy（i0）-（h*h/6）*M0（i0））*（XX（i0+1）-x2）/h+...
（yy（i0+1）-（h*h/6）*M0（i0+1））*（x2-XX（i0））/h;%利用循环求不同区间段的三样条插值多项式
end
Sx%得到不同区间段的三样条插值多项式，用矩阵表示
w1=[XX
（1）:
0.01:
XX
（2）];%将第一段等分曲线段等分成20段，得各等分点横坐标
P1=zeros（1,21）;
fort1=1:
21
x2=w1（t1）;
subs（Sx
（1））;
P1（1,t1）=subs（Sx
（1））;%求得第一段等分曲线段等分成20段的各等分点横坐标对应的纵坐标
end
plot（w1,P1）%绘制利用三样条插值多项式所得的第一等分段曲线
holdon
w2=[XX
（2）:
0.01:
XX（3）];%将第2段等分曲线段等分成20段，得各等分点横坐标
P2=zeros（1,21）;
fort2=1:
21
x2=w2（t2）;
subs（Sx
（2））;
P2（1,t2）=subs（Sx
（2））;%求得第2段等分曲线段等分成20段的各等分点横坐标对应的纵坐标
end
plot（w2,P2）%绘制利用三样条插值多项式所得的第2等分段曲线
holdon
w3=[XX（3）:
0.01:
XX（4）];%将第3段等分曲线段等分成20段，得各等分点横坐标
P3=zeros（1,21）;
fort3=1:
21
x2=w3（t3）;
subs（Sx（3））;
P3（1,t3）=subs（Sx（3））;%求得第3段等分曲线段等分成20段的各等分点横坐标对应的纵坐标
end
plot（w3,P3）%绘制利用三样条插值多项式所得的第3等分段曲线
holdon
w4=[XX（4）:
0.01:
XX（5）];%将第4段等分曲线段等分成20段，得各等分点横坐标
P4=zeros（1,21）;
fort4=1:
21
x2=w4（t4）;
subs（Sx（4））;
P4（1,t4）=subs（Sx（4））;%求得第4段等分曲线段等分成20段的各等分点横坐标对应的纵坐标
end
plot（w4,P4）%绘制利用三样条插值多项式所得的第4等分段曲线
holdon
w5=[XX（5）:
0.005:
XX（6）];%将第5段等分曲线段等分成40段，得各等分点横坐标
P5=zeros（1,41）;
fort5=1:
41
x2=w5（t5）;
subs（Sx（5））;
P5（1,t5）=subs（Sx（5））;%求得第5段等分曲线段等分成40段的各等分点横坐标对应的纵坐标
end
plot（w5,P5）%绘制利用三样条插值多项式所得的第5等分段曲线
holdon
%以下依次绘制其余各等分段的曲线
w6=[XX（6）:
0.005:
XX（7）];
P6=zeros（1,41）;
fort6=1:
41
x2=w6（t6）;
subs（Sx（6））;
P6（1,t6）=subs（Sx（6））;
end
plot（w6,P6）%利用三样条插值多项式绘制第6等分段曲线
holdon
w7=[XX（7）:
0.01:
XX（8）];
P7=zeros（1,21）;
fort7=1:
21
x2=w7（t7）;
subs（Sx（7））;
P7（1,t7）=subs（Sx（7））;
end
plot（w7,P7）%利用三样条插值多项式绘制第7等分段曲线
holdon
w8=[XX（8）:
0.01:
XX（9）];
P8=zeros（1,21）;
fort8=1:
21
x2=w8（t8）;
subs（Sx（8））;
P8（1,t8）=subs（Sx（8））;
end
plot（w8,P8）%利用三样条插值多项式绘制第8等分段曲线
holdon
w9=[XX（9）:
0.01:
XX（10）];
P9=zeros（1,21）;
fort9=1:
21
x2=w9（t9）;
subs（Sx（9））;
P9（1,t9）=subs（Sx（9））;
end
plot（w9,P9）%利用三样条插值多项式绘制第9等分段曲线
holdon
w10=[XX（10）:
0.01:
XX（11）];
P10=zeros（1,21）;
fort10=1:
21
x2=w10（t10）;
subs（Sx（10））;
P10（1,t10）=subs（Sx（10））;
end
plot（w10,P10）%利用三样条插值多项式绘制第10等分段曲线
4.龙格库塔法
%龙格库塔法求3阶常微分方程的初值问题的通用程序
h=input（'请输入步长值h='）%通过键盘输入步长h
a=input（'请输入给定的x的初始点值a='）%通过键盘输入起始点的x值
b=input（'请输入要求近似值的点的x值b='）%通过键盘输入所求近似值点的x值
n=（b-a）/h%求得步长个数
s=input（'请输入m阶微分方程转化为一阶微分方程的后的表达式f（x,y1,y2,y3）=','s'）%输入3阶常微分方程转化为一阶常微分方程后的表达式
f=inline（s）
y=zeros（3,n+1）;%定义一个矩阵来存放循环求得的y1,y2,y3值
k1=zeros（3,n）;%定义一个矩阵来存放循环求得的k1值
k2=zeros（3,n）;%定义一个矩阵来存放循环求得的k2值
k3=zeros（3,n）;%定义一个矩阵来存放循环求得的k3值
k4=zeros（3,n）;%定义一个矩阵来存放循环求得的k4值
A=input（'请将微分方程在初始点a的不同阶的的值按阶数由低到高输入到矩阵中A中构成一个3x1的列矩阵A='）%将给定的初值输入到列矩阵A中
y（:
1）=A;%将y1,y2,y3的初值赋值给矩阵y的第一列
fori=2:
n+1%引入循环
k1（1,i-1）=h*y（2,i-1）;%写出k1矩阵第一行的数学表达式
k1（2,i-1）=h*y（3,i-1）;%写出k1矩阵第二行的数学表达式
k1（3,i-1）=h*f（a,y（1,i-1）,y（2,i-1）,y（3,i-1））;%写出k1矩阵第三行的数学表达式
k2（1,i-1）=h*（y（2,i-1）+k1（2,i-1）/2）;%写出k2矩阵第一行的数学表达式
k2（2,i-1）=h*（y（3,i-1）+k1（3,i-1）/2）;%写出k2矩阵第二行的数学表达式
k2（3,i-1）=h*f（a+（（i-2）*h+h/2）,y（1,i-1）+k1（1,i-1）/2,y（2,i-1）+k1（2,i-1）/2,y（3,i-1）+k1（3,i-1）/2）;%写出k2矩阵第三行的数学表达式
k3（1,i-1）=h*（y（2,i-1）+k2（2,i-1）/2）;%写出k3矩阵第一行的数学表达式
k3（2,i-1）=h*（y（3,i-1）+k2（3,i-1）/2）;%写出k3矩阵第二行的数学表达式
k3（3,i-1）=h*f（a+（（i-2）*h+h/2）,y（1,i-1）+k2（1,i-1）/2,y（2,i-1）+k2（2,i-1）/2,y（3,i-1）+k2（3,i-1）/2）;%写出k3矩阵第三行的数学表达式
k4（1,i-1）=h*（y（2,i-1）+k3（2,i-1））;%写出k4矩阵第一行的数学表达式
k4（2,i-1）=h*（y（3,i-1）+k3（3,i-1））;%写出k4矩阵第二行的数学表达式
k4（3,i-1）=h*f（a+（（i-2）*h+h）,y（1,i-1）+k3（1,i-1）,y（2,i-1）+k3（2,i-1）,y（3,i-1）+k3（3,i-1））;%写出k4矩阵第三行的数学表达式
y（1,i）=y（1,i-1）+（1/6）*（k1（1,i-1）+2*k2（1,i-1）+2*k3（1,i-1）+k4（1,i-1））;%利用循环求得矩阵y的第一行
y（2,i）=y（2,i-1）+（1/6）*（k1（2,i-1）+2*k2（2,i-1）+2*k3（2,i-1）+k4（2,i-1））;%利用循环求得矩阵y的第二行
y（3,i）=y（3,i-1）+（1/6）*（k1（3,i-1）+2*k2（3,i-1）+2*k3（3,i-1）+k4（3,i-1））;%利用循环求得矩阵y的第三行
end
y%得到矩阵y
jg=y（1,n+1）%y（b）的近似值=y（1,n+1）,所求的近似解与解析解值很接近
升级会员 