收藏
下载资源下载资源 加入VIP,免费下载

Jacobi迭代法GaussSeidel迭代法.docx

上传人：b****7 文档编号：8732275 上传时间：2023-02-01 格式：DOCX 页数：18 大小：125.69KB
下载 相关 举报
Jacobi迭代法GaussSeidel迭代法.docx_第1页
第1页 / 共18页
Jacobi迭代法GaussSeidel迭代法.docx_第2页
第2页 / 共18页
Jacobi迭代法GaussSeidel迭代法.docx_第3页
第3页 / 共18页
Jacobi迭代法GaussSeidel迭代法.docx_第4页
第4页 / 共18页
Jacobi迭代法GaussSeidel迭代法.docx_第5页
第5页 / 共18页
点击查看更多>>
下载资源
资源描述

Jacobi迭代法GaussSeidel迭代法.docx

《Jacobi迭代法GaussSeidel迭代法.docx》由会员分享，可在线阅读，更多相关《Jacobi迭代法GaussSeidel迭代法.docx（18页珍藏版）》请在冰豆网上搜索。

Jacobi迭代法GaussSeidel迭代法.docx

Jacobi迭代法GaussSeidel迭代法

Matlab线性方程组的迭代解法（Jacobi迭代法Gauss-Seidel迭代法）实验报告

2008年11月09日星期日12:

49

1.熟悉Jacobi迭代法，并编写Matlab程序matlab程序

按照算法（Jacobi迭代法）编写Matlab程序（Jacobi.m）

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

用Jacobi迭代法求方程组

的解。

输入：

A=[430;33-1;0-14]；

b=[24;30;-24]；

[x,k,index]=Jacobi（A,b,1e-5,100）

输出：

x=

-2.9998

11.9987

-3.0001

k=

100

index=

0

2.熟悉Gauss-Seidel迭代法，并编写Matlab程序

function[v,sN,vChain]=gaussSeidel（A,b,x0,errorBound,maxSp）

%Gauss-Seidel迭代法求解线性方程组

%A-系数矩阵b-右端向量x0-初始迭代点errorBound-近似精度maxSp-最大迭代次数

%v-近似解sN-迭代次数vChain-迭代过程的所有值

step=0;

error=inf;

s=size（A）;

D=zeros（s

（1））;

vChain=zeros（15,3）;%最多能记录15次迭代次数

k=1;

fx0=x0;

fori=1:

s

（1）

D（i,i）=A（i,i）;

end;

L=-tril（A,-1）;

U=-triu（A,1）;

whileerror>=errorBound&step

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

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置：首页 > 初中教育

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1