计算方法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),所求的近似解与解析解值很接近