break;
end
end
X
Gauss-Seidel迭代法运行的解
X=
-0.3941
0.5058
0.7612
-0.7030
3、分析:
Jacobi迭代法称简单迭代法,程序编译比较简单,但是在迭代过程中,对已经算出来的信息未加充分利用,但是该方法算出来的值比较精确,Gauss-Seidel迭代法占用内存小,收敛快,但是精确度稍差。
二、非线性方程的迭代解法
1.用Newton迭代法、割线法求方程
在
附近的正.计算停止的条件为:
;
解:
首先使用使用Newton迭代法
1、MATLAB的编译程序如下:
function[output_args]=Untitled1(input_args)
%UNTITLED1Summaryofthisfunctiongoeshere
%Detailedexplanationgoeshere
%f(x)=x*exp(x)-1
%f'(x)=exp(x)+x*exp(x)
%φ'(x)=x-f(x)/f'(x)
%迭代初值为x=0.5
clc;
clear;
i=1;
Xk=0.5;
Xk1=Xk-(Xk*exp(Xk)-1)/(exp(Xk)+Xk*exp(Xk));
whileabs(Xk1-Xk)>=10^(-6)
Xk=Xk1;
Xk1=Xk-(Xk*exp(Xk)-1)/(exp(Xk)+Xk*exp(Xk));
i=i+1;
end
x=Xk1%x的值即为所求
i
Newton迭代法运行结果
x=
0.5671
i=
4
再运用割线法割线法
2、MATLAB的编译程序如下
function[output_args]=Untitled2(input_args)
%UNTITLED2Summaryofthisfunctiongoeshere
%Detailedexplanationgoeshere
%迭代初值x1=0.5,x2=0.4
%迭代公式为Xk+1=Xk-f(Xk)/((f(Xk)-f(Xk-1))/(Xk-Xk-1))
clc;
clear;
i=1;
Xk1=0.5;
Xk2=0.4;
Xk3=Xk2-(Xk2*exp(Xk2)-1)*(Xk2-Xk1)/((Xk2*exp(Xk2)-1)-(Xk1*exp(Xk1)-1));
whileabs(Xk3-Xk2)>=10^(-6)
Xk1=Xk2;
Xk2=Xk3;
Xk3=Xk2-(Xk2*exp(Xk2)-1)*(Xk2-Xk1)/((Xk2*exp(Xk2)-1)-(Xk1*exp(Xk1)-1));
i=i+1;
end
x=Xk3
i
割线法
运行结果
x=
0.5671
i=
5
3、比较分析:
Newton迭代法是二阶收敛,割线法是在前者基础上变化来的,收敛阶1.618,在该题中前者的迭代步数为4,后者为5,可见Newton迭代法收敛比较快。
三、数值积分
用Romberg方法求
的近似值,要求误差不超过
.
解:
Romberg的编译程序:
function[output_args]=Untitled1(input_args)
%UNTITLED1Summaryofthisfunctiongoeshere
%Detailedexplanationgoeshere
%求解2/sqre(pi)*|exp(-x^2)(01)的值。
clc;
clear;
e=10^(-5);%误差
i=1;%循环控制变量
H0(i)=(f(0)+f
(1))/2;
H1(i)=H0(i)/2+f(0.5)/2;
H1(i+1)=H1(i)*4/3-H0(i)/3;
%两层循环
whileabs(H1(i+1)-H0(i))>=e
Q=0;
a=0;
whilea<=1
a=a+1/2^(i+1);
Q=Q+f(a);
a=a+1/2^(i+1);
end
Q=Q/2^(i+1);
H2
(1)=H1
(1)/2+Q;
forj=2:
i+2
H2(j)=4^(j-1)/(4^(j-1)-1)*H2(j-1)-1/(4^(j-1)-1)*H1(j-1);
end
H0=H1;
H1=H2;
i=i+1;
end
H1(i+1)%存放算法求得的解
function[y]=f(x)
y=2/sqrt(pi)*exp(-x^2);
Romberg方法的程序计算结果
ans=
0.8427
四、插值与逼近
2.已知函数值
0
1
2
3
4
5
6
7
8
9
10
0
0.79
1.53
2.19
2.71
3.03
3.27
2.89
3.06
3.19
3.29
和边界条件:
. 求三次样条插值函数
并画出其图形.
解:
编译程序:
:
在工作空间输入:
clear
closeall
x=0:
10;
y=[00.791.532.192.713.033.272.893.063.193.29];
%求解mi程序
gk=[];
fork=1:
9
gk(k)=3*(y(k+2)-y(k))/2;
end
m0=0.8;
m10=0.2;
A=[20.50000000
0.520.5000000
00.520.500000
000.520.50000
0000.520.5000
00000.520.500
000000.520.50
0000000.520.5
00000000.52];
b1=gk
(1)-0.5*m0;
b9=gk(9)-0.5*m10;
b=[b1gk
(2)gk(3)gk(4)gk(5)gk(6)gk(7)gk(8)b9]';
m=A\b
%图形程序
pp=csape(x,y,'complete',[0.8,0.2]);%第一类边界条件调用函数
xi=0:
0.25:
10;
yi=ppval(pp,xi);
plot(x,y,'o',xi,yi),gridon
title('三次样条插值图(第一类边界条件)')
xlabel('x')
ylabel('y')
结果:
数据结果:
m=
0.77148596314996
0.70405614740014
0.61228944724946
0.38678606360200
0.36056629834254
-0.14905125697216
-0.18436127045388
0.25649633878770
0.05837591530307
图形如下图:
五、常微分方程数值解法
用Euler法与改进的Euler法求解
取步长
计算,画出数值解的图形并与精确解
相比较。
解:
编译程序
function[output_args]=Untitled1(input_args)
%UNTITLED1Summaryofthisfunctiongoeshere
%Detailedexplanationgoeshere
clc;
clear;
h=0.1;
a=0;
b=1;
ya=1;
N=(b-a)/h;
T=linspace(a,b,N+1);
Y=zeros(1,N+1);
Y
(1)=ya;
%Euler法
forj=1:
N
Y(j+1)=Y(j)+h*(Y(j)-T(j)*Y(j)^2);
end
plot(T,Y,'g');gridon;
holdon;
%Euler改进法
forj=1:
N
Y(j+1)=Y(j)+(h/2)*(Y(j)-T(j)*Y(j)^2+(Y(j)+h*(Y(j)-T(j)*Y(j)^2))-T(j+1)*(Y(j)+h*(Y(j)-T(j)*Y(j)^2))^2);
end
plot(T,Y,'r');
holdon;
%精确解
Y=1./(T-1+2*exp(-T));
plot(T,Y);
xlabel('x');
ylabel('f(x)');
title('Euler(绿色)、Euler改进法(红色)及精确解(蓝色)的比较');
分析:
Euler法收敛阶是1阶,改进的Euler算法精度比较高更加逼近真实值。