x0=inv(D)*(L+U)*x0+inv(D)*b;
vChain(k,:
)=x0';
k=k+1;
error=norm(x0-fx0);
fx0=x0;
step=step+1;
end
v=x0;
sN=step;
用Gauss-Seidel迭代法求解上题的线性方程组,取。
输入:
A=[430;33-1;0-14];
b=[24;30;-24];
x0=[0;0;0];
[v,sN,vChain]=gaussSeidel(A,b,x0,0.00001,11)
输出:
v=
0.6169
11.1962
-4.2056
sN=
11
vChain=
6.000010.0000-6.0000
-1.50002.0000-3.5000
4.500010.3333-5.5000
-1.75003.6667-3.4167
3.250010.6111-5.0833
-1.95835.0556-3.3472
2.208310.8426-4.7361
-2.13196.2130-3.2894
1.340311.0355-4.4468
-2.27667.1775-3.2411
0.616911.1962-4.2056
000
000
000
000
s
数值实验
数值实验要求:
数值实验报告内容:
要包含题目、算法公式、完整的程序、正确的数值结果和图形以及相应的误差分析。
在本课程网站上提交数值实验报告的电子文档。
一、为了逼近飞行中的野鸭的顶部轮廓曲线,已经沿着这条曲线选择了一组点。
见下表。
1.对这些数据构造三次自然样条插值函数,并画出得到的三次自然样条插值曲线;
2.对这些数据构造Lagrang插值多项式,并画出得到的Lagrang插值多项式曲线。
x0.91.31.92.12.63.03.94.44.75.06.0
f(x)1.31.51.852.12.62.72.42.152.052.12.25
x7.08.09.210.511.311.612.012.613.013.3
f(x)2.32.251.951.40.90.70.60.50.40.25
1.使用三次样条插值函数csape()求解。
解:
输入:
>>x=[0.91.31.92.12.63.03.94.44.75.06.07.08.09.210.511.311.612.012.613.013.3];
>>y=[1.31.51.852.12.62.72.42.152.052.12.252.32.251.951.40.90.70.60.50.40.25];
>>pp=csape(x,y,'variational',[00])
得到:
pp=
form:
'pp'
breaks:
[0.90001.30001.90002.10002.600033.90004.40004.700056789.200010.500011.300011.60001212.60001313.3000]
coefs:
[20x4double]
pieces:
20
order:
4
dim:
1
再输入:
>>pp.coefs
得到:
ans=
-0.247600.53961.3000
0.9469-0.29720.42081.5000
-2.95641.40731.08681.8500
-0.4466-0.36661.29492.1000
0.4451-1.03650.59342.6000
0.1742-0.5025-0.02222.7000
0.0781-0.0322-0.50342.4000
1.31420.0849-0.47712.1500
-1.58121.2676-0.07132.0500
0.0431-0.15550.26232.1000
-0.0047-0.02610.08082.2500
-0.0244-0.04010.01462.3000
0.0175-0.1135-0.13902.2500
-0.0127-0.0506-0.33581.9500
-0.0203-0.1002-0.53181.4000
1.2134-0.1490-0.73120.9000
-0.83930.9431-0.49290.7000
0.0364-0.0640-0.14130.6000
-0.44800.0014-0.17890.5000
0.5957-0.5361-0.39280.4000
因此,三次样条函数S(x)为
最后输入:
>>xx=0:
0.01:
14;yy=interp1(x,y,xx,'csape');
>>plot(xx,yy);xlabel('x');ylabel('y');
得到图形:
2.Lagrange插值算法:
1)输入N个节点数组;
2)定义初始值和;
3)利用公式
求N次插值基函数;
4)利用多项式求插值函数;
解:
输入:
>>x=[0.9,1.3,1.9,2.1,2.6,3.0,3.9,4.4,4.7,5.0,6.0,7.0,8.0,9.2,10.5,11.3,11.6,12.0,12.6,13.0,13.3];
y=[1.3,1.5,1.85,2.1,2.6,2.7,2.4,2.15,2.05,2.1,2.25,2.3,2.25,1.95,1.4,0.9,0.7,0.6,0.5,0.4,0.25];
>>a=polyfit(x,y,length(x)-1);
>>p=vpa(poly2sym(a),5)
得到数值结果:
p=.30732e-10*x^20+.42770e-8*x^19-.27708e-6*x^18+.11098e-4*x^17-.30784e-3*x^16+.62785e-2*x^15-.97558e-1*x^14+1.1810*x^13-11.297*x^12+86.085*x^11-524.68*x^10+2558.0*x^9-9942.3*x^8+30586.*x^7-73622.*x^6+.13627e6*x^5-.18907e6*x^4+.18913e6*x^3-.12803e6*x^2+52170.*x-9593.4
再输入:
>>y1=polyval(a,x);
plot(x,y1);xlabel('x');ylabel('y');
得到图形:
结果分析:
由以上两图形可以看出三次样条插值的图形较lagrange插值的图形要光滑的多。
因为样条函数插值不仅具有较好的收敛性和稳定性,而且其光滑性也较高,因此样条函数成为了重要的插值工具,其中应用较多的是三次样条插值。
二、对于问题
将h=0.025的Euler法,h=0.05的改进的Euler法和h=0.1的4阶经典的Runge-Kutta法在这些方法的公共节点0.1,0.2,0.3,0.4和0.5处进行比较。
精确解为:
。
1.Euler法
在x,y平面上微分方程的解在曲线上一点(x,y)的切线斜率等于函数的值。
该曲线的顶点设为p,再推进到p(),显然两个顶点p,p的坐标有以下关系
Matlab程序:
function[x,y]=Euler(ydot,x0,y0,h,N)
%ydot为一阶微分方程的函数
%x0,y0为初始条件
%h为区间步长
%N为区间个数
%x为Xn构成的向量,y为Yn构成的向量
x=zeros(1,N+1);y=zeros(1,N+1);
x
(1)=x0;y
(1)=y0;
forn=1:
N
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(ydot,x(n),y(n));
end
解:
取h=0.025,N=20,输入:
>>ydot=inline('y-x^2+1','x','y');
[t,u]=Euler(ydot,0,0.5,0.025,20)
得到数值结果:
t=
Columns1through17
00.02500.05000.07500.10000.12500.15000.17500.20000.22500.25000.27500.30000.32500.35000.37500.4000
Columns18through21
0.42500.45000.47500.5000
u=
Columns1through17
0.50000.53750.57590.61530.65550.69660.73870.78160.82530.87000.91550.96181.00891.05691.10571.15531.2056
Columns18through21
1.25681.30871.36131.4147
即采用Euler法得到:
u(0.1)=0.6555,u(0.2)=0.8253,u(0.3)=1.0089,u(0.4)=1.2056,u(0.5)=1.4147
2.改进Euler方法
改进Euler公式.
y=yn+h/2(f()+f(xn+1,+h*f()))即迭代公式为:
Matlab程序:
function[x,y]=Euler_ad(ydot,x0,y0,h,N)
%改进Euler公式
%ydot为一阶微分方程的函数
%x0,y0为初始条件
%h为区间步长
%N为区间个数
%x为Xn构成的向量,y为Yn构成的向量
x=zeros(1,N+1);y=zeros(1,N+1);
x
(1)=x0;y
(1)=y0;
forn=1:
N
x(n+1)=x(n)+h;
ybar=y(n)+h*feval(ydot,x(n),y(n));
y(n+1)=y(n)+h/2*(feval(ydot,x(n),y(n))+feval(ydot,x(n+1),ybar));
end
解:
取h=0.05,N=10,输入:
>>ydot=inline('y-x^2+1','x','y');
[t,u]=Euler_ad(ydot,0,0.5,0.05,10)
得到数值结果:
t=
00.05000.10000.15000.20000.25000.30000.35000.40000.45000.5000
u=
0.50000.57680.65730.74140.82910.92021.01471.11261.21361.31781.4250
即采用改进的Euler法得到:
u(0.1)=0.6573,u(0.2)=0.8291,u(0.3)=1.0147,u(0.4)=1.2136,u(0.5)=1.4250
3.四阶Runge_Kutta法
求解公式为:
Matlab程序:
function[x,y]=Runge_Kutta4(ydot,x0,y0,h,N)
%标准四阶Runge_Kutta方法
%ydot为一阶微分方程的函数
%x0,y0为初始条件
%h为区间步长
%N为区间个数
%x为Xn构成的向量,y为Yn构成的向量
x=zeros(1,N+1);y=zeros(1,N+1);
x
(1)=x0;y
(1)=y0;
forn=1:
N
x(n+1)=x(n)+h;
k1=h*feval(ydot,x(n),y(n));
k2=h*feval(ydot,x(n)+1/2*h,y(n)+1/2*k1);
k3=h*feval(ydot,x(n)+1/2*h,y(n)+1/2*k2);
k4=h*feval(ydot,x(n)+h,y(n)+k3);
y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);
end
解:
取h=0.1,N=5,输入:
>>ydot=inline('y-x^2+1','x','y');
[t,u]=Runge_Kutta4(ydot,0,0.5,0.1,5)
得到数值结果:
t=
00.10000.20000.30000.40000.5000
u=
0.50000.65740.82931.01511.21411.4256
结果比较:
tEuler法改进Euler法四阶Runge_Kutta精确解
0.10.65550.65730.65740.6574
0.20.82530.82910.82930.8293
0.31.00891.01471.01511.0151
0.41.20561.21361.21411.2141
0.51.41471.42501.42561.4256
由以上结果可以看出改进的Euler法较Euler法计算精度有所提高,但还不是十分精确。
四阶Runge_Kutta法具有非常高的精度,事实上,在求解微分方程初值问题,四阶Runge_Kutta法是单步长中最优秀的方法,通常都用该方法求解实际问题。
三、
用Newton迭代法求方程的根时,分别取初始值,;
用Newton迭代法求方程时,分别取初始值,;
算法:
(1)取初始点x0最大迭代次数N和精度要求ε,k=0.
(2)如果f’(xk)=0,则停止计算;否则计算
Xk+1=xk-f(xk)/f’(xk)
(3)若|xk+1-xk|<ε,则停止计算。
(4)若k=N,则停止计算;否则置k=k+1,转
(2)。
Matlab程序:
function[x_star,index,it]=Newton(fun,x,ep,it_max)
%求解非线性方程的Newton法
%fun(x)为需要求根的函数,第一个分量是函数值,第二个分量是导数值
%fun=inline('[x^3-x-1,3*x^2-1]')当f(x)=x^3-x-1;
%x为初始点
%ep为精度,缺省值为1e-5
%it_max为最大迭代次数,缺省值100
%x_star为当迭代成功时输出方程的根,失败时输出最后的迭代值
%index为指标变量,index=1表示迭代成功index=0表示失败
%it为迭代次数
ifnargin<4it_max=100;end
ifnargin<3ep=1e-5;end
index=0;k=1;
whilek<=it_max
x1=x;f=feval(fun,x);
ifabs(f
(2))x=x-f
(1)/f
(2);
ifabs(x-x1)k=k+1;
end
x_star=x;it=k;
解:
(1)由于f(x)=arctan(x),f’(x)=1/1+x2,取初始值x0=1.45,输入
>>fun=inline('[atan(x),1/(1+x^2)]');
>>[x_star,index,it]=Newton(fun,1.45)
得到数值结果:
x_star=1.6281e+004
index=0
it=7
取初始值x0=0.5,输入
>>fun=inline('[atan(x),1/(1+x^2)]');
>>[x_star,index,it]=Newton(fun,0.5)
得到数值结果:
x_star=0
index=1
it=4
输入x=-3:
0.05:
3;y=atan(x);
plot(x,y);xlabel('x');ylabel('y');
得到图形:
(2)由于f(x)=x3-x-3=0,f’(x)=3x2-1,取初始值x0=0.0,输入
>>fun=inline('[x^3-x-3,3*x^2-1]');
>>[x_star,index,it]=Newton(fun,0.0)
得到数值结果:
x_star=-0.0074
index=0
it=101
取初始值x0=2.0,输入
>>fun=inline('[x^3-x-3,3*x^2-1]');
>>[x_star,index,it]=Newton(fun,2.0)
得到数值结果:
x_star=1.6717
index=1
it=4
输入x=0:
0.05:
3;y=x.^3-x-3;
plot(x,y);xlabel('x');ylabel('y');
得到图形:
结果分析:
从以上结果可以看出初值的选取与Newton法的收敛很有关系。
初值选的不好,Nexton法将不收敛,它的收敛性是在跟a附近讨论的。
选取初值时一定要十分小心。
四、用Jacobi迭代和SOR法分别求解线性方程组(教科书第86页算例2),并验证、输出SOR法的松弛因子w和对应的迭代次数。
1.Jacobi迭代法
Jacobi迭代法的算法为:
Matlab程序:
function[x,k,index]=Jacobi(A,b,ep,it_max)
%求解线性方程组的Jacobi迭代法
%A为系数矩阵
%b为方程组右端项
%ep为精度要求,缺省值1e-5
%it_max为最大迭代次数,缺省值100
%x为方程组的解
%k为迭代次数
%index为指标变量index=1表示迭代收敛到指定要求index=0表示迭代失败。
ifnargin<4it_max=100;end
ifnargin<3ep=1e-5;end
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);index=1;
while1
fori=1:
n
y(i)=b(i);
forj=1:
n
ifj~=i
y(i)=y(i)-A(i,j)*x(j);
end
end
ifabs(A(i,i))<1e-10|k==it_max
index=0;return;
end
y(i)=y(i)/A(i,i);
end
ifnorm(y-x,inf)break;
end
x=y;k=k+1;
end
function[A,b]=shuru
%自己定义输入矩阵A和向量b的函数
n=15;
fori=1:
n
ifi==1|i==15z(i)=3;
elsez(i)=2;
end
forj=1:
n
ifj==iA(i,j)=4;
elseifj==i+1|j==i-1A(i,j)=-1;
elseA(i,j)=0;
end
end
end
b=z';
解:
输入:
>>[A,b]=shuru;ep=1e-6;
>>[x,k,index]=Jacob