当前位置：文档库 › 计算电磁学之FDTD算法的MATLAB语言实现

计算电磁学之FDTD算法的MATLAB语言实现

South China Normal University

课程设计实验报告

课程名称：计算电磁学

指导老师：

专业班级：2014级电路与系统

姓名：

学号：

FDTD算法的MATLAB语言实现

摘要：时域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接对麦克斯韦方

程作差分处理,用来解决电磁脉冲在电磁介质中传播和反射问题的算法。其基本思想是:FDTD计算域空间节点采用Yee元胞的方法,同时电场和磁场节点空间与时间上都采用交错抽样；把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区,这两个区域是以连接边界相连接,最外边是采用特殊的吸收边界,同时在这两个边界之间有个输出边界,用于近、远场转换;在连接边界上采用连接边界条件加入入射波,从而使得入射波限制在总场区域；在吸收边界上采用吸收边界条件,尽量消除反射波在吸收边界上的非物理性反射波。

本文主要结合FDTD算法边界条件特点,在特定的参数设置下，用MATLAB语言进行编程，在二维自由空间TEz网格中，实现脉冲平面波。

关键词：FDTD；MATLAB；算法

1绪论

1.1课程设计背景与意义

20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(TLM)和时域积分方程法等。其中FDTD是一种已经获得广泛应用并且有很大发展前景的时域数值计算方法。时域有限差分（FDTD）方法于1966年由K.S.Yee提出并迅速发展，且获得广泛应用。K.S.Yee用后来被称作Yee氏网格的空间离散方式，把含时间变量的Maxwell旋度方程转化为差分方程，并成功地模拟了电磁脉冲与理想导体作用的时域响应。但是由于当时理论的不成熟和计算机软硬件条件的限制，该方法并未得到相应的发展。20世纪80年代中期以后，随着上述两个条件限制的逐步解除，FDTD便凭借其特有的优势得以迅速发展。它能方便、精确地预测实际工程中的大量复杂电磁问题，应用范围几乎涉及所有电磁领域，成为电磁工程界和理论界研究的一个热点。目前，FDTD日趋成熟，并成为分析大部分实际电磁问题的首选方法。

1.2时域有限差分法的发展与应用

经过四十多年的发展，FDTD已发展成为一种成熟的数值计算方法。在发展

过程中，几乎都是围绕几个重要问题展开的，即数值稳定性、计算精度、数值色散、激励源技术以及开域电磁问题的吸收边界条件等。数值稳定和计算精度对任何一种数值计算方法都是至关重要的。其中激励源的设计和引入也是FDTD的一个重要任务。目前，应用最广泛的激励源引入技术是总场/散射场体系。对于散射问题，通常在FDTD计算空间中引入连接边界，它将整个计算空间划分为内部的总场区和外部的散射场区，如图1.1。利用Huygens原理，可以在连接边界处引入入射场，使入射场的加入变得简单易行。

图1.1总场/散射场区

2FDTD法基本原理

2.1Maxwell方程和Yee氏算法

根据电磁场基本方程组的微分形式，若在无源空间，其空间中的煤质是各向同性、线性和均匀的，即煤质的参数不随时间变化且各向同性，则Maxwell旋度

方程可写成：

式中，E是电场强度，单位为伏/米（V/m）；H是磁场强度，单位为安/米（A/m）；ε表示介质介电系数，单位为法拉/米（F/m）；μ表示磁导系数，单位为亨利/米（H/m）；σ表示介质电导率，单位为西门子/米（S/m）；σ

表示导磁率，单位

为欧姆/米（Ω/m）。

K.S.Yee在1966年建立了如图2.1所示的空间网格，这就是著名的Yee氏元胞网格。

图2.1Yee氏网格及其电磁场分量分布

电场和磁场被交叉放置,电场分量位于网格单元每条棱的中心,磁场分量位于网格单元每个面的中心,每个磁场(电场)分量都有4个电场(磁场)分量环绕。这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足,而且还允许旋度方程在空间上进行中心差分运算,同时也满足了法拉第电磁感应定律和安培环路积分定律,也可以很恰当地模拟电磁波的实际传播过程。

2.2时域有限差分法相关技术

2.2.1数值稳定性问题

FDTD方程是一种显式差分方程，在执行时，存在一个重要的问题：即算法的稳定性问题。这种不稳定性表现为在解显式方程时，随着时间步数的继续增加，计算结果也将无限制地增加。Taflove于1975年对Yee氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。数值解是否稳定主要取决于时间步长Δt与空间步长Δx、Δy、Δz的关系。对于非均匀媒质构成的计算空间选用如下的稳定性条件：

上式是空间和时间离散之间应当满足的关系，又称为Courant稳定性条件。若采用均匀立方体网格：Δx=Δy=Δz=Δs

其中，v为计算空间中的电磁波的最大速度。

2.2.2数值色散

FDTD方程组是对Maxwell旋度方程进行差分近似，在进行数值计算时，将会在计算网格中引起数字波模的色散，即在FDTD网格中，电磁波的相速与频率有关，电磁波的相速度随波长、传播方向及变量离散化的情况不同而改变。这种关系由非物理因素引起，且色散将导致非物理因素引起的脉冲波形畸变、人为的各向异性和虚假折射等现象。为获得理想的色散关系，问题空间分割应按照小于正常网格的原则进行。一般选取的最大空间步长为Δmax=λmin/20，λmin为所研究范围内电磁波的最小波长。由上分析说明，数值色散在用FDTD法分析电磁场传播中的影响是不可能避免的，但我们可以尽可能的减小数值色散的影响。

2.2.3离散网格的确定

无论是简单目标还是复杂目标，在进行FDTD离散时网格尺寸的确定，除了受计算资源的限制不可能取得很小外，还需要考虑以下几个因素：1．目标离散精确度的要求。网格应当足够小以便能精确模拟目标几何形状和电磁参数。

2．FDTD方法本身的要求。主要是考虑色散误差的影响。设网格为立方体Δ

x=Δy=Δz=δ，所关心频段的频率上限为f

max ，对应波长为f

min

，则考虑FDTD

的数值色散要求

δ≤λ

min

通常N≥10。上式是根据已知所关心频率上限情况下来确定FDTD网格尺寸的；反之，若给定，则FDTD计算结果可用的上限频率也随之确定。

2.3吸收边界条件

由时域有限差分法的基本原理可知，在利用时域有限差分法研究电磁场时，需在全部问题空间建立Yee氏网格空间，并存储每个单元网格上任一时间步的六个场分量用于下一时间步的计算。而在对于辐射、散射这类开放系统的实际研究中，不可能有无限大的存储空间。因此，必须在某处将网格空间截断，且在截断边界网格点处运用特殊的场分量计算方法，使得向边界面行进的波在边界处保持“外向行进”特性、无明显的反射现象，并且不会使内部空间的场产生畸变，从而用有限网格空间模拟电磁波在无界空间中传播的情况。具有这种功能的边界条件称之为吸收边界条件，或辐射边界条件，或网格截断条件，如图2.2所示

图2.2附加截断边界使计算区域变为有限域

从FDTD的基本差分方程组可以看出，在截断边界面上切向场分量的计算需要利用计算空间以外的电磁场分量，因此FDTD基本差分方程对这些截断边界面上的场分量失效。如何处理截断边界上的场分量，使之与需要考虑的无限空间有尽量小的差异，是FDTD中必须很好解决的一个重要问题。实际上，这是要求在误差可容忍的范围内，计算空间中的外向波能够顺利通过截断边界面而不引起波的明显反射，使有限计算空间的数值模拟与实际情况趋于一致，对外向波而言，就像在无限大空间中传播一样。所以，需要一种截断边界网格处的特殊计算方法，它不仅要保证边界场计算的必要精度，而且还要大大消除非物理因素引起的波反射，使得用有限的网格空间就能模拟电磁波在无限空间中的传播。但是如果处理不当，截断边界面可能造成较大反射，构成数值模拟误差的一部分，甚至可能造成算法不稳定。

加于截断边界场分量符合上述要求的算法就称为吸收边界条件（Absorbing Boundary Conditions）。

2.4FDTD计算所需时间步的估计

为了使计算达到稳定,通常计算所需要时间步按照电磁波往返穿越FDTD计算区对角线3～5次来估计。若FDTD计算区总元胞数为N，则对角线上元胞约为（三维)。按照Courant稳定条件，设计算中心区Δt=δ/（2c），即穿

越对角线一次需要时间步为。总计算时间步约为需。对于二维

情况则约为。或者说，计算时间步大约等于FDTD计算区对角线上元胞数目的12～20倍。实际上，计算所需时间步还与散射体具体形状、结构有关。

图2.3给出了应用FDTD分析电磁场问题时的程序流程图

图2.3FDTD程序流程图

3课程设计要求及MATLAB语言实现

3.1课程设计要求

采用总场/散射场技术，在二维自由空间TEz网格中，实现脉冲平面波。总场区域100*100个网格。

周边散射场区域20个网格。

散射场区域外采用PML吸收边界条件。

平面波0度角入射。使用高斯脉冲，中心频率为1GHz,带宽为1/e。

采用正方形Yee网格，离散步长Δ=1.5cm，时间离散步长为Δ=Δ/2c.

3.2MATLAB程序

%%%%%%%%此程序实现自由空间的平面波%%%%%%%% %%%%%%%%此程序采用总场/散射场技术%%%%%%%% %%%%%%%%此程序使用PML设置吸收边界条件%%%%%%%%

KE=128;%真空区域网格数

IE=KE;

JE=KE;

T=0;

NSTEP=300;%迭代次数

e=2.71828;

%%%%%%%%初始化%%%%%%%%

dz=zeros(KE,KE);

ez=zeros(KE,KE);

hx=zeros(KE,KE);

hy=zeros(KE,KE);

ihx=zeros(KE,KE);

ihy=zeros(KE,KE);

ga=ones(KE,KE);

ez_inc=zeros(1,JE);

hx_inc=zeros(1,JE);

%%%%%%%%PML参数%%%%%%%%

gi2=ones(1,IE);gi3=gi2;

fi1=zeros(1,IE);

fi2=ones(1,IE);fi3=fi2;

gj2=ones(1,JE);gj3=gj2;

fj1=zeros(1,JE);

fj2=ones(1,JE);fj3=fj2;

for i=1:8%设置PML层中的参数

xnum=npml-i;

xd=npml;

xxn=xnum/xd;

xn=0.33*xxn^3;

gi2(i)=1.0/(1.0+xn);

gi2(IE-1-i)=1.0/(1.0+xn);

gi3(i)=(1.0-xn)/(1.0+xn);

gi3(IE-1-i)=(1.0-xn)/(1.0+xn);

xxn=(xnum-0.5)/xd;

xn=0.25*xxn^3;

fi1(i)=xn;

fi1(IE-2-i)=xn;

fi2(i)=1.0/(1.0+xn);

fi2(IE-2-i)=1.0/(1.0+xn);

fi3(i)=(1.0-xn)/(1.0+xn);

fi3(IE-2-i)=(1.0-xn)/(1.0+xn); end

for j=1:8

xnum=npml-j;

xd=npml;

xxn=xnum/xd;

xn=0.33*xxn^3;

gj2(j)=1.0/(1.0+xn);

gj2(JE-1-j)=1.0/(1.0+xn);

gj3(j)=(1.0-xn)/(1.0+xn);

gj3(JE-1-j)=(1.0-xn)/(1.0+xn);

xxn=(xnum-0.5)/xd;

xn=0.25*xxn^3;

fj1(j)=xn;

fj1(JE-2-j)=xn;

fj2(j)=1.0/(1.0+xn);

fj2(JE-2-j)=1.0/(1.0+xn);

fj3(j)=(1-xn)/(1.0+xn);

fj3(JE-2-j)=(1.0-xn)/(1.0+xn); end

%%%%%%%%总场/散射场边界%%%%%%%%

ib=IE-ia-1;

ja=28;

jb=JE-ja-1;

ez_inc_low_m1=0;

ez_inc_low_m2=0;

ez_inc_high_m1=0;

ez_inc_high_m2=0;

%%%%%%%%信号源%%%%%%%%

t0=80;%脉冲高度

spread=12;%脉冲宽度

%%%%%%%%网格%%%%%%%%

ddx=0.015;%%%%%%%%离散步长%%%%%%%%

dt=ddx/(2*3e8);%%%%%%%%计算时间离散步长%%%%%%%%

epsz=8.851*1e-12;%%%%%%%%自由空间的介电常数%%%%%%%%

%%%%%%迭代求解电场与磁场%%%%%%%%

for n=1:NSTEP

a=1;

T=T+1;

for j=2:JE

ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));

end

%%%%%%%%入射波缓冲区%%%%%%%%

ez_inc(1)=ez_inc_low_m2;

ez_inc_low_m2=ez_inc_low_m1;

ez_inc_low_m1=ez_inc(2);

ez_inc(JE-1)=ez_inc_high_m2;

ez_inc_high_m2=ez_inc_high_m1;

ez_inc_high_m1=ez_inc(JE-2);

%%%%%%%%计算dx区域%%%%%%%%

for i=2:KE

for j=2:KE

dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx (i,j)+hx(i,j-1));

end

%%%%%%%%输入高斯脉冲信号源%%%%%%%%

pulse=exp(-0.5*(t0-T)^2/spread^2);%%%%%%%%

ez_inc(3)=pulse;

%%%%%%%%入射波DZ参数%%%%%%%%

for i=ia:ib

dz(i,ja)=dz(i,ja)+0.5*hx_inc(ja-1);

dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb);

end

%%%%%%%%计算电场dx区域%%%%%%%%

for i=1:KE

for j=1:KE

ez(i,j)=ga(i,j)*dz(i,j);

end

%%%%%%%%设置电场EZ边缘为0,作为PML的参数%%%%%%%%

for j=1:JE-1

ez(1,j)=0;

ez(IE,j)=0;

end

for i=1:IE-1

ez(i,1)=0;

ez(i,JE)=0;

end

for j=1:JE-1

hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1));

end

%%%%%%%%计算磁场hx区域%%%%%%%%

for i=1:KE

for j=1:KE-1

curl_e=ez(i,j)-ez(i,j+1);

ihx(i,j)=ihx(i,j)+fi1(i)*curl_e;

hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j));

end

%%%%%%%%入射磁场Hx区域参数%%%%%%%%

for i=ia:ib

hx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja);

hx(i,jb)=hx(i,jb)-0.5*ez_inc(jb);

end

%%%%%%%%计算磁场hy区域%%%%%%%%

for i=1:KE-1

for j=1:KE

curl_e=ez(i+1,j)-ez(i,j);

ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;

hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j));

end

%%%%%%%%入射磁场hy区域参数%%%%%%%%

for j=ja:jb

hy(ia-1,j)=hy(ia-1,j)-0.5*ez_inc(j);

hy(ib,j)=hy(ib,j)+0.5*ez_inc(j);

end

%%%%%%%%图形展示%%%%%%%%

m=moviein(500);

x=1:KE;

y=1:KE;

[X,Y]=meshgrid(x,y);

surf(X,Y,ez);

axis([0KE0KE-11]);

set(gcf,'Color','white','Number','off','Name',

sprintf('Simulation FDTD2D,Iteration=%i',n));

title(sprintf('t=%.3f nsec',n*dt*1e9))

xlabel('x');

ylabel('y');

zlabel('Ez');

m(a)=getframe(gcf);

a=a+1;

end

下图为该程序实现的效果图：

4.结语

以上结合FDTD和MATLAB在自由空间中实现了平面波，所编MATLAB程序简洁明了，运行效率也较高。FDTD法在电磁场数值分析方面有很大的优越性，而MATLAB具有强大的数据处理和图形处理功能，可以快速地编出高效高质量的程序。将二者的优势有效地结合起来，可以将算法迅速程序化，并获得很好的数据处理结果。但由于我的编程思路不是很规范，所以在实现过程中，遇到了很多的问题，比如在完全匹配层吸收效果上稍有差距，但已经基本上达到了预期的效果。

参考文献：

[1]葛德彪,闫玉波.电磁场时域有限差分方法[M].西安:西安电子科技大学出版社,2005.

[2]MATLAB模拟的电磁学时域有限差分法.[美]Atef Elsherbeni Veysel Demir著.

喻志远译.国防工业出版社，2012.

[3]Engquist B,Majda A.Absorbing Boundary Conditions for the Numerical Simulation of Waves. Mathematics of Computation.1977.

[4](美)Duane Hanselman,Bruce Littlefield著.精通Matlab7[M].朱仁峰,译.北京:清华大学出版社,2006.5.

[5]数值解法和MATLAB实现与应用.Gerald Recktenward著.伍卫国万群张辉等译.

[6]MATLAB数字计算与仿真.周品编著.电子工业出版社2013.

数值计算方法实验指导(Matlab版)

《数值计算方法》实验指导（Matlab 版）肇庆学院数学与统计学学院计算方法课程组

1. 实验名称实验1 算法设计原则验证(之相近数相减、大数吃小数和简化计算步骤) 2. 实验题目有效数字的损失. 123 )与1000个较小的数(3 10 15)的和，验证大数吃小数的现象. (3)分别用直接法和秦九韶算法计算多项式 P(x) a 0x n a 1x n 1 在x =1.00037 处的值?验证简化计算步骤能减少运算时间. n 1 对于第(3)题中的多项式P (x )，直接逐项计算需要n (n 1) 2 1 次乘法和n 次加法，使用秦九韶算法 P(x) (((a °x ajx a 2)x a . 则只需要n 次乘法和n 次加法. 3. 实验目的验证数值算法需遵循的若干规则. 4. 基础理论设计数值算法时，应避免两个相近的数相减、防止大数吃小数、简化计算步骤减少运算次数以减少运算时间并降低舍入误差的积累. 两相近的数相减会损失有效数字的个数，用一《数值计算方法》实验 1报告班级： 20xx 级 XXXXx 班学号： 20xx2409xxxx 姓名： XXX 成绩: ⑴取 z 1016，计算z 1 Z 和 1/（、z 1 Z)，验证两个相近的数相减会造成 (2)按不同顺序求一个较大的数( a n 1 X a n

个大数依次加小数，小数会被大数吃掉，乘法运算次数太多会增加运算时间. 5.实验环境操作系统：Win dows xp ；程序设计语言：Matlab 6.实验过程 (1)直接计算并比较； (2)法1 :大数逐个加1000个小数，法2 :先把1000个小数相加再与大数加； (3)将由高次项到低次项的系数保存到数组A[n]中，其中n为多项式次数. 7.结果与分析 (1)计算的~1V Z = _______________________________ ,1/( ~1 < z) ____________________ . 分析： (2)123逐次加1000个3 10 6的和是_________________________ ，先将1000个3 10 6相加，再用这个和与123相加得_______________________ . 分析： (3)计算__________ 次的多项式：直接计算的结果是___________________ ，用时___________________ ; 用秦九韶算法计算的结果是____________________ ，用时 ___________________ 分析：

科学计算与MATLAB语言考试答案

1 单选(2分) 利用MATLAB求解科学计算问题的优势是（）。得分/总分 ? A. 算法最优 ? B. 不需要编写程序 ? C. 程序执行效率高 ? D. 编程效率高正确答案：D你没选择任何选项 2 单选(2分) 在MATLAB命令行窗口输入命令时，可使用续行符，其写法是（）。得分/总分 ? A. 省略号（…） ? B. 分号（;） ? C. 三个小数点（…） ? D. 百分号（%）正确答案：C你没选择任何选项 3

下列语句执行后，D的值为（）。 1.A=[1:3;4:6]; 2.D=sub2ind(size(A),[1,1],[2,3]) 得分/总分 ? A. 3 6 ? B. 2 5 ? C. 3 5 ? D. 4 5 正确答案：C你没选择任何选项 4 单选(2分) ceil(-2.1)+floor(-2.1)+fix(-2.1)的结果为（）。得分/总分 ? A. -7 ? B. -6 ? C. -5 ? D. -9 正确答案：A你没选择任何选项 5

下列语句执行后，x的值是（）。 1.log=1:5; 2.x=log(1) 得分/总分 ? A. ? B. 1 ? C. 数学常数e ? D. 报错正确答案：B你没选择任何选项 6 单选(2分) 下列语句执行后，c的值是（）。 1.ch=['abcdef';'123456']; 2.c=char(ch(2,4)-1) 得分/总分 ? A. '4' ? B. 4 ? C. '3' ? D. 3

7 单选(2分) 产生和A同样大小的全0矩阵的函数是（）。得分/总分 ? A. zero(size(A)) ? B. zeros(size(A)) ? C. size(zero(A)) ? D. size(zeros(A)) 正确答案：B你没选择任何选项 8 单选(2分) 语句x=speye(5)==eye(5)执行后，则下列说法中正确的是（）。得分/总分 ? A. x是5阶全1矩阵，且采用稀疏存储方式 ? B. x是5阶全1矩阵，且采用完全存储方式 ? C. x是5阶单位矩阵，且采用稀疏存储方式 ? D. x是5阶单位矩阵，且采用完全存储方式

计算方法_全主元消去法_matlab程序

%求四阶线性方程组的MA TLAB程序 clear Ab=[0.001 2 1 5 1; 3 - 4 0.1 -2 2; 2 -1 2 0.01 3; 1.1 6 2.3 9 4];%增广矩阵 num=[1 2 3 4];%未知量x的对应序号 for i=1:3 A=abs(Ab(i:4,i:4));%系数矩阵取绝对值 [r,c]=find(A==max(A(:))); r=r+i-1;%最大值对应行号 c=c+i-1;%最大值对应列号 q=Ab(r,:),Ab(r,:)=Ab(i,:),Ab(i,:)=q;%行变换 w=Ab(:,c),Ab(:,c)=Ab(:,i),Ab(:,i)=w;%列变换 n=num(i),num(i)=num(c),num(c)=n;%列变换引起未知量x次序变化for j=i:3 Ab(j+1,:)=-Ab(j+1,i)*Ab(i,:)/Ab(i,i)+Ab(j+1,:);%消去过程 end end %最后得到系数矩阵为上三角矩阵 %回代算法求解上三角形方程组 x(4)=Ab(4,5)/Ab(4,4); x(3)=(Ab(3,5)-Ab(3,4)*x(4))/Ab(3,3); x(2)=(Ab(2,5)-Ab(2,3)*x(3)-Ab(2,4)*x(4))/Ab(2,2); x(1)=(Ab(1,5)-Ab(1,2)*x(2)-Ab(1,3)*x(3)-Ab(1,4)*x(4))/Ab(1,1); for s=1:4 fprintf('未知量x%g =%g\n',num(s),x(s)) end %验证如下 %A=[0.001 2 1 5 1; 3 -4 0.1 -2 2;2 -1 2 0.01 3; 1.1 6 2.3 9 4]; %b=[1 2 3 4]'; %x=A\b; %x1= 1.0308 %x2= 0.3144 %x3= 0.6267 %x4= -0.0513

数值分析Matlab作业

数值分析编程作业

2012年12月第二章 14.考虑梯形电阻电路的设计，电路如下：电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组： 12 123 234 345 456 567 678 78 22/ 2520 2520 2520 2520 2520 2520 250 i i V R i i i i i i i i i i i i i i i i i i i i -= -+-= -+-= -+-= -+-= -+-= -+-= -+= 这是一个三对角方程组。设V=220V，R=27Ω，运用追赶法，求各段电路的电流量。Matlab程序如下： function chase () %追赶法求梯形电路中各段的电流量 a=input('请输入下主对角线向量a='); b=input('请输入主对角线向量b='); c=input('请输入上主对角线向量c='); d=input('请输入右端向量d='); n=input('请输入系数矩阵维数n='); u(1)=b(1); for i=2:n l(i)=a(i)/u(i-1); u(i)=b(i)-c(i-1)*l(i); end y(1)=d(1); for i=2:n y(i)=d(i)-l(i)*y(i-1); end x(n)=y(n)/u(n); i=n-1; while i>0 x(i)=(y(i)-c(i)*x(i+1))/u(i); i=i-1; end x 输入如下:

请输入下主对角线向量a=[0,-2,-2,-2,-2,-2,-2,-2]; 请输入主对角线向量b=[2,5,5,5,5,5,5,5]; 请输入上主对角线向量c=[-2,-2,-2,-2,-2,-2,-2,0]; 请输入方程组右端向量d=[220/27,0,0,0,0,0,0,0]; 请输入系数矩阵阶数n=8 运行结果如下： x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477 第三章 14.试分别用(1)Jacobi 迭代法；(2)Gauss-Seidel 迭代法解线性方程组 1234510123412191232721735143231211743511512x x x x x ?????? ??????---????????????=--?????? --?????? ??????---?????? 迭代初始向量 (0)(0,0,0,0,0)T x =。（1）雅可比迭代法程序如下： function jacobi() %Jacobi 迭代法 a=input('请输入系数矩阵a='); b=input('请输入右端向量b='); x0=input('请输入初始向量x0='); n=input('请输入系数矩阵阶数n='); er=input('请输入允许误差er='); N=input('请输入最大迭代次数N='); for i=1:n for j=1:n if i==j d(i,j)=a(i,j); else d(i,j)=0; end end end m=eye(5)-d\a; %迭代矩阵 g=d\b; x=m*x0+g; k=1; while k<=N %进行迭代 for i=1:5 if max(abs(x(i)-x0(i))) >er x=m*x+g; k=k+1;

蚁群算法TSP问题matlab源代码

function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta ,Rho,Q) %%===================================================== ==================== %% ACATSP.m %% Ant Colony Algorithm for Traveling Salesman Problem %% ChengAihua,PLA Information Engineering University,ZhengZhou,China %% Email:aihuacheng@https://www.wendangku.net/doc/9e5354461.html, %% All rights reserved %%------------------------------------------------------------------------- %% 主要符号说明 %% C n个城市的坐标，n×4的矩阵 %% NC_max 最大迭代次数 %% m 蚂蚁个数 %% Alpha 表征信息素重要程度的参数 %% Beta 表征启发式因子重要程度的参数 %% Rho 信息素蒸发系数 %% Q 信息素增加强度系数 %% R_best 各代最佳路线 %% L_best 各代最佳路线的长度 %%===================================================== ==================== %%第一步：变量初始化 n=size(C,1);%n表示问题的规模（城市个数） D=zeros(n,n);%D表示完全图的赋权邻接矩阵 for i=1:n for j=1:n if i~=j D(i,j)=max( ((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5,min(abs(C(i,3)-C(j,3)),144- abs(C(i,3)-C(j,3))) );%计算城市间距离 else D(i,j)=eps; end D(j,i)=D(i,j); end end Eta=1./D;%Eta为启发因子，这里设为距离的倒数 Tau=ones(n,n);%Tau为信息素矩阵 Tabu=zeros(m,n);%存储并记录路径的生成 NC=1;%迭代计数器 R_best=zeros(NC_max,n);%各代最佳路线

数值分析的matlab实现

第2章牛顿插值法实现参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55. 求牛顿插值多项式和差商的MA TLAB 主程序： function[A,C,L,wcgs,Cw]=newpoly(X,Y) n=length(X);A=zeros(n,n);A(:,1) =Y'; s=0.0;p=1.0;q=1.0;c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1; end C=A(n,n);b=poly(X(n));q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C);Q=poly2sym(q1); syms M wcgs=M*Q/c1;Cw=q1/c1; （1）保存名为newpoly.m 的M 文件（2）输入MA TLAB 程序 >> X=[242,243,249,250]; >> Y=[13.681,13.526,13.098,13.095]; >> [A,C,L,wcgs,Cw]=newpoly(X,Y) 输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其 ) ()()1(ξ+n n f x R 的系数向量Cw 。 A = 13.6810 0 0 0 13.5260 -0.1550 0 0 13.0980 -0.0713 0.0120 0 13.0950 -0.0030 0.0098 -0.0003 C = 1.0e+003 *

中南大学材料学院科学计算与MATLAB考试题库

练习题 1.求函数在指定点的数值导数 x=sym('x'); >> y=[x x.^2 x.^3;1 2*x 3*x.^2;0 2 6*x]; >> x=1; >> eval(diff(y)) ans = 1 2 3 0 2 6 0 0 6 >> x=2; >> eval(diff(y)) ans = 1 4 12 0 2 12 0 0 6 >> x=3; >> eval(diff(y)) ans = 1 6 27 0 2 18 0 0 6 2.求下列函数导数（1） x=sym('x'); >> y=x^10+10^x+(log(10))/log(x); >> diff(y) ans = 10*x^9+10^x*log(10)-2592480341699211/1125899906842624/log(x)^2/x （2） x=sym('x');

>> y=log(1+x); >> x=1; >> eval(diff(y,2)) %在x=1的条件下对y表达式求两次导数后导函数的值 ans = -0.2500 3.用数值方法求下列积分首先先讲一下trapz的用法，如下题 t=0:0.001:1; >> y=t; >> trapz(t,y) ans = 0.5000 （1） >> x=1:0.01:5; >> y=(x.^2).*sqrt(2*x.^2+3); >> trapz(x,y) ans = 232.8066 (2) x=pi/4:0.01:pi/3; >> y=x./(sin(x).^2); >> trapz(x,y) ans = 0.3810 第三题拟合曲线题 x=[0:0.1:1]; >> y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; >> a=polyfit(x,y,2); >> x=[0.05:0.2:1.05]; >> y=a(3)+a(2)*x+a(1)*x.^2 %注意x要在y前先赋值，不然y不会运行为最新的x对呀的y值 y =

第2讲 matlab的数值分析

第二讲MATLAB的数值分析 2-1矩阵运算与数组运算矩阵运算和数组运算是MATLAB数值运算的两大类型，矩阵运算是按矩阵的运算规则进行的，而数组运算则是按数组元素逐一进行的。因此，在进行某些运算（如乘、除）时，矩阵运算和数组运算有着较大的差别。在MATLAB中，可以对矩阵进行数组运算，这时是把矩阵视为数组，运算按数组的运算规则。也可以对数组进行矩阵运算，这时是把数组视为矩阵，运算按矩阵的运算规则进行。 1、矩阵加减与数组加减矩阵加减与数组加减运算效果一致，运算符也相同，可分为两种情况：（1）若参与运算的两矩阵（数组）的维数相同，则加减运算的结果是将两矩阵的对应元素进行加减，如 A=[1 1 1;2 2 2;3 3 3]; B=A; A+B ans= 2 2 2 4 4 4 6 6 6 （2）若参与运算的两矩阵之一为标量（1*1的矩阵），则加减运算的结果是将矩阵（数组）的每一元素与该标量逐一相加减，如 A=[1 1 1;2 2 2;3 3 3]; A+2 ans= 3 3 3 4 4 4 5 5 5 2、矩阵乘与数组乘（1）矩阵乘矩阵乘与数组乘有着较大差别，运算结果也完全不同。矩阵乘的运算符为“*”，运算是按矩阵的乘法规则进行，即参与乘运算的两矩阵的内维必须相同。设A、B为参与乘运算的 =A m×k B k×n。因此，参与运两矩阵，C为A和B的矩阵乘的结果，则它们必须满足关系C m ×n 算的两矩阵的顺序不能任意调换，因为A*B和B*A计算结果很可能是完全不一样的。如：A=[1 1 1;2 2 2;3 3 3]; B=A;

A*B ans= 6 6 6 12 12 12 18 18 18 F=ones(1,3); G=ones(3,1); F*G ans 3 G*F ans= 1 1 1 1 1 1 1 1 1 （2）数组乘数组乘的运算符为“.*”，运算符中的点号不能遗漏，也不能随意加空格符。参加数组乘运算的两数组的大小必须相等（即同维数组）。数组乘的结果是将两同维数组（矩阵）的对应元素逐一相乘，因此，A.*B和B.*A的计算结果是完全相同的，如： A=[1 1 1 1 1;2 2 2 2 2;3 3 3 3 3]; B=A; A.*B ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 B.*A ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 由于矩阵运算和数组运算的差异，能进行数组乘运算的两矩阵，不一定能进行矩阵乘运算。如 A=ones(1,3); B=A; A.*B ans= 1 1 1 A*A ???Error using= =>

数值分析的MATLAB程序

列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1);

end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x')

王能超计算方法——算法设计及MATLAB实现课后代码

第一章插值方法 1.1Lagrange插值 1.2逐步插值 1.3分段三次Hermite插值 1.4分段三次样条插值第二章数值积分 2.1 Simpson公式 2.2 变步长梯形法 2.3 Romberg加速算法 2.4 三点Gauss公式第三章常微分方程德差分方法 3.1 改进的Euler方法 3.2 四阶Runge-Kutta方法 3.3 二阶Adams预报校正系统 3.4 改进的四阶Adams预报校正系统第四章方程求根 4.1 二分法 4.2 开方法 4.3 Newton下山法 4.4 快速弦截法第五章线性方程组的迭代法 5.1 Jacobi迭代 5.2 Gauss-Seidel迭代 5.3 超松弛迭代 5.4 对称超松弛迭代第六章线性方程组的直接法 6.1 追赶法 6.2 Cholesky方法 6.3 矩阵分解方法 6.4 Gauss列主元消去法

第一章插值方法 1.1Lagrange插值计算Lagrange插值多项式在x=x0处的值. MATLAB文件：（文件名：Lagrange_eval.m）function [y0,N]= Lagrange_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Lagrange插值多项式在x0处的值 %N是Lagrange插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; for j=1:m if j~=i; N(i)=N(i)*(x0-X(j))/(X(i)-X(j)); end end y0=y0+Y(i)*N(i); end 用法》X=[…];Y=[…]; 》x0= ; 》[y0,N]= Lagrange_eval(X,Y,x0) 1.2逐步插值计算逐步插值多项式在x=x0处的值. MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Neville逐步插值多项式在x0处的值 m=length(X); P=zeros(m,1); P1=zeros(m,1); P=Y; for i=1:m P1=P; k=1; for j=i+1:m k=k+1;

科学计算与matlab1.5

单元测验已完成成绩：分 1 【单选题】 MATLAB一词来自（）的缩写。 A、 Mathematica Laboratory B、 Matrix Laboratory C、 MathWorks Lab D、 Matrices Lab 我的答案：B得分：分 2 【单选题】下列选项中能反应MATLAB特点的是（）。 A、算法最优 B、不需要写程序 C、程序执行效率高 D、编程效率高

我的答案：D得分：分单元测验已完成成绩：分 1 【单选题】当在命令行窗口执行命令时，如果不想立即在命令行窗口中输出结果，可以在命令后加上（）。 A、冒号（:） B、逗号（,） C、分号（;） D、百分号（%）我的答案：C得分：分 2 【单选题】 fix(264/100)+mod(264,10)*10的值是（）。 A、 86 B、 62 C、 423 D、

42 我的答案：D得分：分 3 【单选题】在命令行窗口输入下列命令后，x的值是（）。 >> clear >> x=i*j A、不确定 B、 -1 C、 1 D、 i*j 我的答案：B得分：分 4 【单选题】使用语句x=linspace(0,pi,6)生成的是（）个元素的向量。 A、 8 B、 7 C、 6

D、 5 我的答案：C得分：分 5 【单选题】 ceil的结果为（）。 A、 -2 B、 -3 C、 1 D、 2 我的答案：A得分：分 6 【单选题】 eval('sqrt(4)+2')的值是（）。 A、 sqrt(4)+2 B、 4 C、 2 D、

2+2 我的答案：B得分：分 7 【单选题】已知a为3×5矩阵，则执行完a(:,[2,4])=[]后（）。 A、 a变成行向量 B、 a变为3行2列 C、 a变为3行3列 D、 a变为2行3列我的答案：C得分：分 8 【单选题】在命令行窗口输入以下命令 >> A=[1:3;4:6]; >> D=sub2ind(size(A),[1,1],[2,3]) D的值为（）。 A、 3 6 B、 2 5 C、 4 5

(整理)matlab16常用计算方法.

常用计算方法 1．超越方程的求解一超越方程为 x (2ln x – 3) -100 = 0 求超越方程的解。 [算法]方法一：用迭代算法。将方程改为 01002ln()3 x x =- 其中x 0是一个初始值，由此计算终值x 。取最大误差为e = 10-4，当| x - x 0| > e 时，就用x 的值换成x 0的值，重新进行计算；否则| x - x 0| < e 为止。 [程序]P1_1abs.m 如下。 %超越方程的迭代算法 clear %清除变量 x0=30; %初始值 xx=[]; %空向量 while 1 %无限循环 x=100/(2*log(x0)-3); %迭代运算 xx=[xx,x]; %连接结果 if length(xx)>1000,break ,end %如果项数太多则退出循环（暗示发散） if abs(x0-x)<1e-4,break ,end %当精度足够高时退出循环 x0=x; %替换初值 end %结束循环 figure %创建图形窗口 plot(xx,'.-','LineWidth',2,'MarkerSize',12)%画迭代线'.-'表示每个点用.来表示，再用线连接 grid on %加网格 fs=16; %字体大小 title('超越方程的迭代折线','fontsize',fs)%标题 xlabel('\itn','fontsize',fs) %x 标签 ylabel('\itx','fontsize',fs) %y 标签 text(length(xx),xx(end),num2str(xx(end)),'fontsize',fs)%显示结果 [图示]用下标作为自变量画迭代的折线。如P0_20_1图所示，当最大误差为10-4时，需要迭代19次才能达到精度，超越方程的解为27.539。 [算法]方法二：用求零函数和求解函数。将方程改为函数 100()2ln()3f x x x =-- MATLAB 求零函数为fzero ，fzero 函数的格式之一是 x = fzero(f,x0) 其中，f 表示求解的函数文件，x0是估计值。fzero 函数的格式之二是 x = fzero(f,[x1,x2])

同济大学数值分析matlab编程题汇编

MATLAB 编程题库 1.下面的数据表近似地满足函数2 1cx b ax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像. 625 .0718.0801.0823.0802.0687.0606.0356.0995 .0628.0544.0008.0213.0362.0586.0931.0i i y x ---- 解： x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y; a=z(1); b=z(2); c=z(3); xh=-1:0.1:1; yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*')

2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数 10 的近似根，并写出调用方式：精度为10 解： >> edit gexianfa.m function [x iter]=gexianfa(f,x0,x1,tol) iter=0; while(norm(x1-x0)>tol) iter=iter+1; x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end >> edit f.m function v=f(x) v=x.*log(x)-1; >> edit g.m function z=g(y) z=y.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8

科学计算与MATLAB语言(第四课)

第四讲绘图功能

作为一个功能强大的工具软件，Matlab 具有很强的图形处理功能，提供了大量的二维、三维图形函数。由于系统采用面向对象的技术和丰富的矩阵运算，所以在图形处理方面即常方便又高效。

4.1 二维图形一、plot函数函数格式：plot(x,y)其中x和y为坐标向量函数功能：以向量x、y为轴，绘制曲线。【例1】在区间0≤X≤2 内，绘制正弦曲线Y=SIN（X），其程序为： x=0:pi/100:2*pi; y=sin(x); plot(x,y)

一、plot函数【例2】同时绘制正、余弦两条曲线Y1=SIN（X）和Y2=COS（X），其程序为： x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,x,y2) plot函数还可以为plot(x,y1,x,y2，x,y3，…)形式，其功能是以公共向量x为X轴，分别以y1，y2，y3，…为Y轴，在同一幅图内绘制出多条曲线。

一、plot函数（一）线型与颜色格式：plot(x,y1,’cs’,...) 其中c表示颜色，s表示线型。【例3】用不同线型和颜色重新绘制例4.2图形，其程序为：x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,'go',x,y2,'b-.') 其中参数'go'和'b-.'表示图形的颜色和线型。g表示绿色，o表示图形线型为圆圈；b表示蓝色，-.表示图形线型为点划线。

一、plot函数（二）图形标记在绘制图形的同时，可以对图形加上一些说明，如图形名称、图形某一部分的含义、坐标说明等，将这些操作称为添加图形标记。 title(‘加图形标题'); xlabel('加X轴标记'); ylabel('加Y轴标记'); text(X,Y,'添加文本');

matlab用于计算方法的源程序

1、Newdon迭代法求解非线性方程 function [x k t]=NewdonToEquation(f,df,x0,eps) %牛顿迭代法解线性方程 %[x k t]=NewdonToEquation(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:原函数，定义为内联函数 ?:函数的倒数，定义为内联函数 %x0:初始值 %eps:误差限 % %应用举例： %f=inline('x^3+4*x^2-10'); ?=inline('3*x^2+8*x'); %x=NewdonToEquation(f,df,1,0.5e-6) %[x k]=NewdonToEquation(f,df,1,0.5e-6) %[x k t]=NewdonToEquation(f,df,1,0.5e-6) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1) if nargin==3 eps="0".5e-6; end tic; k=0; while 1 x="x0-f"(x0)./df(x0); k="k"+1; if abs(x-x0) < eps || k >30 break; end x0=x; end t=toc; if k >= 30 disp('迭代次数太多。'); x="0"; t="0"; end

2、Newdon迭代法求解非线性方程组 function y="NewdonF"(x) %牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量 y(1,1)=x(1).^2-10*x(1)+x(2).^2+8; y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8; return; function y="NewdonDF"(x) %牛顿迭代法解非线性方程组的测试函数的导数 y(1,1)=2*x(1)-10; y(1,2)=2*x(2); y(2,1)=x(2).^+1; y(2,2)=2*x(1).*x(2)-10; return; 以上两个函数仅供下面程序的测试 function [x k t]=NewdonToEquations(f,df,x0,eps) %牛顿迭代法解非线性方程组 %[x k t]=NewdonToEquations(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:方程组（事先定义） ?:方程组的导数（事先定义） %x0:初始值 %eps:误差限 % %说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示 % %应用举例： %x0=[0,0];eps=0.5e-6; %x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

《科学计算与MATLAB》期末大作业

杭州电子科技大学信息工程学院《科学计算与MATLAB》期末大作业

给出程序、图、作业分析，程序需加注释。 1. 试编写名为fun.m 的MATLAB 函数，用以计算下述的值： ?? ? ??-<->=t t n t t t n t f 的)4/sin()(si 对所有)4/sin(其他情况)sin(的)4/sin()(si 对所有)4/sin()(ππππ 绘制t 关于函数f(t)的图形，其中t 的取值范围为ππ66≤≤-t ，间距为10/π。 function y=fun()%定义函数 % t=-6*pi:pi/10:6*pi; %定义变量范围 y = (sin(pi/4)).*(sin(t)>sin(pi/4))+(sin(-pi/4)).*(sin(t)=sin(-pi/4)));%函数表示 plot(t,y); %画图 end

2.解以下线性方程组 ??? ??=+=++=--3 530 42231 321321x x x x x x x x A=[2 -1 -1;1 1 4;3 0 5];%输入矩阵 B=[2;0;3]; %输入矩阵 X = A\B %计算结果 3.已知矩阵? ? ??? ???? ???=44434241 3433323124232221 14131211A 求： (1)A(2:3,2:3) (2)A(:,1:2) (3)A(2:3,[1,3]) (4)[A,[ones(2,2)；eye(2)]]

A=[11 12 13 14;21 22 23 24;31 32 33 34;41 42 43 44];%输入矩阵A(2:3,2:3) %输出矩阵 A(:,1:2) %输出矩阵 A(2:3,[1,3]) %输出矩阵 [A,[ones(2,2);eye(2)]] %输出矩阵

数值计算方法与Matlab样卷答案

腹有诗书气自华《数值计算方法与Matlab 》样卷答案一．填空题：（每空3分,共42分） 1． 8，6105.0-? 。 2．)(3)1(2)1(1)(3)1(2)1(1)(3)1(3)(3)(2)1(1)(3)(2)1(1)(2)1(2)(3)(2)(1)(3)(2)(1)(1)1(1)1(22)22()1()1(222)1()222(k k k k k k k k k k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x x x x x x x x x x ωωωωωωωωωω ωωωω-+--=---?+=+--+-=---?+=++--=+--?+=+++++++++， )2,1(∈ω。 3．],[1b a C S m -∈。4. 1e 2e ---x ，???==-=?--? ,3,2,1,0;0,e 1d )(e 110k k x x g k x ，正交投影。 5. 2阶，6阶。 6．10.6658，10.9521，10.9501。 7. 4002.2)00.1(=ε，4030.2)01.1(=ε。二．解下列各题：（每题9分，共36分） 1．解：令)1(2 3+=t x ，（2分）则??-+++=+1123 02 dt )1(25.21)1(49d 1t t x x x ???++++???++-+-≈22)6.01(25.21)6.01(9525.219 8)6.01(25.21)6.01(9549 （8分） 210631.10≈ （9分） 2．解：记系数矩阵为A, 对增广矩阵[]b A |作初等行运算, ??????????--401533933112??????????--==5.55.115 .35.405.75.401125.1,5.11,31,2l l ??????????---=45.114005.75.4011212,3l , 所以13-=x ,2)5.75.1(5.4112=-=x x ,1)1(2 1321=-+-=x x x ,即方程组的解为 [1,2,-1]T . （4分）故系数矩阵A 的LU 分解为???? ??????--???????????---=4005.75.40112115.1015.1001A 。（6分）

0计算方法及MATLAB实现简明讲义课件PPS8-1欧拉龙格法

第8章常微分方程初值问题数值解法 8.1 引言 8.2 欧拉方法 8.3 龙格-库塔方法 8.4 单步法的收敛性与稳定性 8.5 线性多步法

8.1 引言考虑一阶常微分方程的初值问题 00(,),[,],(). y f x y x a b y x y '=∈=（1.1）（1.2）如果存在实数，使得 121212(,)(,).,R f x y f x y L y y y y -≤-?∈（1.3）则称关于满足李普希茨(Lipschitz ）条件，称为的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页)

计算方法及MATLAB 实现所谓数值解法，就是寻求解在一系列离散节点 )(x y <<<<<+121n n x x x x 上的近似值 . ,,,,,121+n n y y y y 相邻两个节点的间距称为步长. n n n x x h -=+1 如不特别说明，总是假定为定数， ),2,1( ==i h h i 这时节点为 . ) ,2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”. 即求解过程顺着节点排列的次序一步一步地向前推进. 00(,),[,], (). y f x y x a b y x y '=∈=

描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算的递推公式. 1+n y 一类是计算时只用到前一点的值，称为单步法. 1+n y n y 另一类是用到前面点的值， 1+n y k 11,,,+--k n n n y y y 称为步法. k 其次，要研究公式的局部截断误差和阶，数值解与精确解的误差估计及收敛性，还有递推公式的计算稳定性等问题. n y )(n x y 首先对方程离散化，建立求数值解的递推公式. ),(y x f y ='