MATLAB语言课程论文基于MATLAB的电磁场数值图像分析

基于MATLAB的电磁场数值分析应用

[摘要] MATLAB使用计算机进行电磁场数值分析已成为电磁场的工程开发、科研和教学的重要手段。编程实现从电磁场微分方程到有限元求解全过程需要很好的理论基础和编程技巧，难度较高。该文介绍了电磁场数值分析的基本理论并通过几个实例介绍了使用MATLAB 实现电磁场偏微分方程的有限元解法。实验结果表明这一方法具有操作简单明了!运算速度快，计算误差可控制等优点

[关键词电磁场数值分析MATLAB 麦克斯韦方程

一、问题的提出

电磁学是物理学的一个分支，是研究电场和电磁的相互作用现象。电磁学从原来互相独立的两门科学(电学、磁学)发展成为物理学中一个完整的分支学科，主要是基于电流的磁效应和变化的磁场的电效应的发现。这两个实验现象，加上麦克斯韦关于变化电场产生磁场的假设，奠定了电磁学的整个理论体系，发展了对现代文明起重大影响的电工和电子技术。

针对电磁场学习理论性强、概念抽象等特点，利用Matlab强大的数值计算和图形技术，通过具体实例进行仿真，绘制相应的图形，使其形象化，便于对其的理解和掌握。将Matlab引入电磁学中，利用其可视化功能对电磁学实验现象进行计算机模拟，可以提高学习效率于学习积极性，使学习效果明显。

通过Matlab软件工具，对点电荷电场、线电荷产生的电位、平面上N个电荷之间的库仑引力、仿真电荷在变化磁场中的运动等问题分别给出了直观形象的的仿真图和数值分析，形实现了可视化学习，丰富了学习内容，提高了对电磁场理论知识的兴趣。从而更好地解决电磁场中数值分析的问题。

二、电磁场数值解法

麦克斯韦方程组是电磁场理论的基础，也是电磁场数值分析的出发点。它的微分形式方程：

ρ=??=????-

=????+=

??→

→→→→→→D B t B

E t D H J c 0

（1）

式中磁场强度电通密度电场强度磁感应强度。电磁场中各种场量之间的关系由

媒质的特性确定。在各向同性媒质中，由下列结构方程组确定

→

→→→→

→===E J H

B E

D σμε (2)

为获得电磁场问题的唯一解!除上述方程组之外尚需给出定解条件，对静态

场和稳态场只需加边界条件，对时变场还需另加初始条件。边界条件包括：

（1）第一类边界条件是给定边界上的值，其中是边界点的函数或常数；

（2）第二类边界条件给定边界上法向导数的值；

（3）第三类边界条件给定边界值与法向导数的线性组合

()()P P n f f 21=+???? (3)

根据麦克斯韦方程组和结构方程组!在静电场中，以电位 为求解对象，在各向

同性介质中，电位满足泊松方程：

ερφ=?

2 (4) 在恒定磁场中，取矢量磁位为求解对象，令有

J c

A μ-=?2 (5) 考虑电磁波动方程： 002

22222=-=-??????t H H t

E E εμεμ (6) 在正弦稳态条件下,上式可分别导出亥姆霍兹方程:

00

2222=+=+???

???H H E E εμεμ?? (7) 如上述几个例子,对于不同的电磁场实际应用问题求解可以得到对应的电磁

场偏微分方程,直接用解析法求解这些方程组往往会遇到很多困难甚至无法求解,电磁场数值分析方法已成为求解电磁场问题的重要方法。

数值分析方法将原来连续的场域离散化求解， 再用离散点的结果近似逼近

连续场域的解&常用的电磁场数值分析方法包括有限差分法、边界元法和有限元

法。 有限差分法是以差分原理为基础的一种数值计算方法，即用差分方程代替

偏微分方程，把要求解的边值问题转化为一组相应的差分问题，将求解区域划

分，求解差分方程组从而得出各网格单元的场值。

这种方法的特点是方法简单，网格划分容易，但对不规则边界处理不便，

网格划分缺乏灵活性。边界元法以麦克斯韦积分方程为基础，它采用分步积分

如格林定理等在一定条件下把该积分方程转化为关于边界的积分方程，并据此

进行离散获得对应的代数方程!求出场域中变量的数值。 它的特点是将数值法

与解析法相结合， 在数学上起到降维的作用，减少了计算量，但对非线性情况

失去了高精度特点，有局限性。

有限元法由于单元定义灵活，处理边界条件容易，具有正定对称系数矩阵

而占据主导地位&有限元法是根据变分原理和离散化而取得近似解的方法。是先

从偏微分方程边值问题出发，找出一个能量泛函的积分式，并令其在满足第一

类边界条件的前提下取极值，即构成条件变分问题& 例如与上式对应的二维泊

松场第三类边界条件下的泛函极值问题为：

()min 21221221'222=??

? ??-+??????????-??????????+=????? ??????? ?????dl dxdy F f f y x L s ?ε?ε??ρ?? f L 02=? (8)

与上式对应的亥姆霍兹方程泛函为：

()(){}dl dxdy F L n

D k ???--=

222221???? (9)

这个条件变分问题是和偏微分方程边值问题等价的。有限元法便是以条件

变分问题为对象来求解电磁场问题。在求解过程中，将场的求解区域剖分成有

限个单元，因此在单元中构造出插值函数，将插值函数代入能量泛函的积分

式，再把泛函离散化成多元函数。通过多元函数极值求极值方法得到一个代数

方程组。最后由此方程组求解得到数值解。下面举几例论证：

例证1

设圆线圈的中心为O ，半径为R ，放置于y-z 平面，线圈通过的电流为I0，

如右图所示。用毕奥－萨伐尔定律计算载流圆线圈在z=0处x-y 平面上的磁场分

布。

解题分析

根据毕奥－萨伐尔定律，

3

d d 4L I r μ?=π?l r B (10)

线圈上任一点处的电流元在x-y 平面上一点P 产生的元磁场为dB 。在编制

程序时，将电流环分为N 段，每一小段视为一电流元，然后求出每一电流元在

观察点处的磁场分量，求出总磁场，最后叠加。

MATLAB 程序如下：

clear all %清除

R=input('请输入圆环半径,R='); %定义输入变量

I0=input('请输入电流,I0='); %定义输出变量

mu0=4*pi*1e-7; C0=mu0/(4*pi); %归并常数

N=20; %电流环分段数

x=linspace(-3,3,N); y=x; %确定观测点范围

theta0=linspace(0,2*pi,N+1); %环的圆周角分段

theta1=theta0(1:N);

y1=R*cos(theta1); z1=R*sin(theta1); %环各段矢量的起始坐标y1,z1

theta2=theta0(2:N+1);

y2=R*cos(theta2); z2=R*sin(theta2); %环各段矢量的终点坐标y2,z2

xc=0; yc=(y2+y1)./2; zc=(z2+z1)./2;

%计算环各段矢量中点的三个坐标分量xc,yc,zc

dlx=0;dly=y2-y1;dlz=z2-z1;

%计算环各段矢量dl 的三个长度分量，其中x1=x2=0。

NGx=N; NGy=NGx; %网格线数

for i=1:NGy %循环计算各网点上的B （x,y ）值

for j=1:NGx

rx=x(j)-xc; ry=y(i)-yc; rz=0-zc;

%计算径矢r 的3个长度分量，r 在z=0的面上。

r3=sqrt(rx.^2+ry.^2+rz.^2).^3; %计算r3

dlXr_x=dly.*rz-dlz.*ry; %计算叉乘dl ×r 的x 和y 分量，z 分量为0

dlXr_y=dlz.*rx-dlx.*rz;

Bx(i,j)=sum(C0*I0.*dlXr_x./r3); %把环各段产生的磁场分量累加

By(i,j)=sum(C0*I0.*dlXr_y./r3);

B=(Bx.^2+By.^2).^0.5; %计算B 的大小

end

end

subplot(1,2,1), quiver(x,y,Bx,By), %画矢量场图

hold on

plot(0,1,'ro',0,-1,'bo'),

xlabel('x'),ylabel('y'), %修饰图形，标注坐标轴

axis([-3,3,-3,3]),

subplot(1,2,2)

mesh(x,y,B);axis([-3,3,-3,3,0,1e-4]) %画磁场大小分布图

xlabel('x'),ylabel('y'),zlabel('B')

运行该程序，例如，在命令窗中的R和I0的提示后分别键入1和100，运行结果如

图1所示

图1 载流圆线圈的磁场分布度及位置曲线

例证2

一对相同的圆形线圈，彼此平行而共轴。设两线圈内的电流都是I，且回

绕方向一致，线圈的半径为R，二者的间距为a（当a=R时，称为亥姆霍兹线

圈），求轴线附近的磁场分布。

解题分析

本题是把观测区域取在两线圈之间的小范围内。线圈B生成的左边的磁场等

于线圈A的左边磁场。因此，A、B两线圈在中间部分的合成磁场等于A线圈的

右磁场与其左磁场平移R后的和。

MATLAB程序如下：

clear all %清除

mu0=4*pi*1e-7;C0=mu0/(4*pi); %归并常数

I0=5.0;R=1;

NGx=21;NGy=21; %设定网格线数

x=linspace(-1,1,NGx); %确定观测点范围

y=linspace(-1,1,NGy);

N=20; %电流环分段数

theta0=linspace(0,2*pi,N+1); %环的圆周角分段

theta1=theta0(1:N);

y1=R*cos(theta1); z1=R*sin(theta1); %环各段矢量的起始坐标y1,z1

theta2=theta0(2:N+1);

y2=R*cos(theta2); z2=R*sin(theta2); %环各段矢量的终点坐标y2,z2

dlx=0;dly=y2-y1;dlz=z2-z1;

%计算环各段矢量dl的三个长度分量，其中x1=x2=0。

xc=0; yc=(y2+y1)/2; zc=(z2+z1)/2;

for i=1:NGy

for j=1:NGx

rx=x(j)-xc; ry=y(i)-yc; rz=0-zc; %计算径矢r的3个长度分量,r在z=0的面上。

r3=sqrt(rx.^2+ry.^2+rz.^2).^3;

dlXr_x=dly.*rz-dlz.*ry; %计算叉乘dl×r的x和y分量，z分量为0

dlXr_y=dly.*rx-dlx.*rz;

Bx(i,j)=sum(C0*I0*dlXr_x./r3); %把环各段产生的磁场分量累加

By(i,j)=sum(C0*I0*dlXr_y./r3);

end

end

Bax=Bx(:,11:21)+Bx(:,1:11); %把x<0区域内的磁场平移，叠加到x>0区域

Bay=By(:,11:21)+By(:,1:11);

subplot(1,2,1),mesh(x(11:21),y,Bax);xlabel('x');ylabel('y');zlabe l('B');

subplot(1,2,2),quiver(x,y,Bx,By,1.5),

axis('square'),axis([-1,1,-1,1]),

xlabel('x');ylabel('y');

运行结果如图2所示。可以看出，在轴线附近磁场大小均匀且沿x方向。

图2 亥姆霍兹线圈轴线附近Bx 在x-y 平面上的分布及矢量场

例证3

设带电粒子质量为m ，带电量为，电场强度E 沿方向，磁感应强度B 沿qyz

方向. 则带电粒子在均匀电磁场中的运动微分方程为

=-=-===???????

?z x

v y y

v x m qB E m q m qB E m q m qB m qB x y (11) 令， 则上面微分方程可化作：

()()()()()()()()()()()06,6524,4342,21==-====dt dy y dt dy y m

qB E m q dt dy y dt dy y m

qB dt dy y dt dy (12)

选择和为参量，就可以分别研究0≠E ，0=B 和，等情况. 编写MATLAB 程序

如下：

clear %清除

syms w x y z t B E m q; %定义变量

E=input('E=');B=input('B='); %输入E 和B 值

[x,y,z]=dsolve('D2x=q*B/m*Dy','D2y=q*E/m-q*B/m*Dx','D2z=0','x(0)

=0','y(0)=0','z(0)=0','Dx(0)=0.01','Dy(0)=6','Dz(0)=0.01') ;

%初始条件取x(0)=y(0)=z(0)=0,Dx(0)=0.01,Dy(0)=6,Dz(0)=0.01 q=1.6e-2; m=0.02; %赋值

X=subs([x y z]); x=X(1),y=X(2),z=X(3),

ezplot3(X(1),X(2),X(3)) %绘图函数调用

运行上述程序，例如，取E=4, B=8可得下列特解并给出图运行结果如图3所示。研究时可以采用不同的初始条件和不同的参量观察不同的现象。运算特解如下：

x =

(t*exp(1))/8 - (100*i*exp(1) - 8*i + 4800)/(10240*exp((32*i*t)/5)) - (exp((32*i*t)/5)*(8*i - 100*i*exp(1) + 4800))/10240 + 15/16

y =

(5*exp(1))/256 + (i*(100*i*exp(1) - 8*i + 4800))/(10240*exp((32*i*t)/5)) - (i*exp((32*i*t)/5)*(8*i - 100*i*exp(1) + 4800))/10240 - 1/640

z =

t/100

图3现有E=4, B=8参数运行结果

例证4

如右图所示，求垂直于无限长载流直导线的平

面内磁感应强度的分布。

解题分析

设场点P 的位置为000x y z ++i j k ，电流元位置为x y z ++i j k ，电流元矢量为

d (d dy d )I I x z =++l i j k 。由此，场点P 相对于电流元的位置矢量为

000()()()r x x y y z z =-+-+-i j k

利用行列式计算I d l ×r ，可写为 000d d d d l r x

y z x x y y z z

?=---i

j k （13）

也可利用MATLAB 中的det 命令函数来求该行列式，MATLAB 程序如下：

syms dx dy dz x0 x y0 y z0 z; %定义变量

dl=[dx,dy,dz]; %定义行列式

r=[x0-x,y0-y, z0-z];

d1cr=cross(dl,r) %求dl ×r 的积

运行结果为

d1cr = [ dy*(z0-z)-dz*(y0-y), dz*(x0-x)-dx*(z0-z),

dx*(y0-y)-dy*(x0-x)]

即[][][]000000d ()d ()d ()d ()d ()d ()d l r z z y y y z x x z z z x y y x x x y ?=---+---+---i j k 又，r 的大小为 222000()()()r x x y y z z =-+-+- 设载流导体通过坐标原点

垂直于 x -y 平面放置，电流元 I d l 沿z 轴正向，场点P 位于x -y 平面上。对本

题目而言，dx=dy=0,x=y=0, z 0=0, 矢量叉乘积为00d d d y z x z ?=-+l r i j ，r 的大小

为 22200r x y z =++

由毕奥－萨伐尔定律 0

3d d 4L I r μ?=π?l r

B （14）

有 002223/200d d 4()x I y z B x y z μ-=

π++； 002223/200d d 4()y I x z B x y z μ=π++ 002223/200d 4()x I

y z B x y z μ∞-∞-=π++?； 002223/200d 4()

y I x z B x y z μ∞-∞=π++? 22x y B B B =+ （15）

MATLAB 程序

1、 用符号运算求B 的表达式

syms C0 I z x y r r0; %定义变量

Bx=C0.*I.*int(-y./(x.^2+y.^2+z.^2).^(3/2),'z',-inf,inf) By=C0.*I.*int(x./(x.^2+y.^2+z.^2).^(3/2),'z',-inf,inf)

B=(Bx.^2+By.^2).^0.5

运行结果：

Bx =-2*C0*I*y/(x^2+y^2)^(3/2)/(1/(x^2+y^2))^(1/2)

By =2*C0*I*x/(x^2+y^2)^(3/2)/(1/(x^2+y^2))^(1/2)

B =(4*C0^2*I^2*y^2/(x^2+y^2)^2+4*C0^2*I^2*x^2/(x^2+y^2)^2)^(1/2) 即 0222

2d 4()x I y z B x y μ=-π+； 02222d 4()y I

x z B x y μ=π+； 022

24I

B x y μ=π+ （16） 2、绘制磁场大小分布图和矢量场图

x=-0.5:0.05:0.5;y=x; %赋值

I=input('请输入电流I='); %设置输入

mu0=4*pi*1e-7; C0=mu0/(4*pi); %设置函数

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

Bx =-2.*C0.*I.*Y./(X.^2+Y.^2).^(3./2.)./(1./(X.^2+Y.^2)).^(1./2); By =2.*C0.*I.*X./(X.^2+Y.^2).^(3./2)./(1./(X.^2+Y.^2)).^(1./2); B=(4.*C0.^2.*I.^2.*Y.^2./(X.^2+Y.^2).^2+4.*C0.^2.*I.^2.*X.^2./(X.^2+Y.^2).^2).^(1./2);

subplot(1,2,1) %选择1×2区域中的1号区

quiver(X,Y,Bx,By,2), axis([-0.5,0.5,-0.5,0.5]), axis('square'), subplot(1,2,2) %选择1×2区域中的2号区

mesh(X,Y,B) %绘图

运行该程序，在命令窗中的提示后键入I0值 (例如，取I0＝100A)，便得到图4所示图形。

MATLAB数据分析与多项式计算(M)

第7章 MATLAB数据分析与多项式计算 6.1 数据统计处理 6.2 数据插值 6.3 曲线拟合 6.4 离散傅立叶变换 6.5 多项式计算 6.1 数据统计处理 6.1.1 最大值和最小值 MATLAB提供的求数据序列的最大值和最小值的函数分别为max 和min，两个函数的调用格式和操作过程类似。 1．求向量的最大值和最小值 求一个向量X的最大值的函数有两种调用格式，分别是： (1) y=max(X)：返回向量X的最大值存入y，如果X中包含复数元素，则按模取最大值。 (2) [y,I]=max(X)：返回向量X的最大值存入y，最大值的序号存入I，如果X中包含复数元素，则按模取最大值。 求向量X的最小值的函数是min(X)，用法和max(X)完全相同。 例6-1 求向量x的最大值。 命令如下： x=[-43,72,9,16,23,47]; y=max(x) %求向量x中的最大值 [y,l]=max(x) %求向量x中的最大值及其该元素的位置 2．求矩阵的最大值和最小值 求矩阵A的最大值的函数有3种调用格式，分别是： (1) max(A)：返回一个行向量，向量的第i个元素是矩阵A的第i 列上的最大值。 (2) [Y,U]=max(A)：返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。 (3) max(A,[],dim)：dim取1或2。dim取1时，该函数和max(A)完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。 求最小值的函数是min，其用法和max完全相同。

例6-2 分别求3×4矩阵x中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。 3．两个向量或矩阵对应元素的比较 函数max和min还能对两个同型的向量或矩阵进行比较，调用格式为： (1) U=max(A,B)：A,B是两个同型的向量或矩阵，结果U是与A,B 同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。 (2) U=max(A,n)：n是一个标量，结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。 min函数的用法和max完全相同。 例6-3 求两个2×3矩阵x, y所有同一位置上的较大元素构成的新矩阵p。 6.1.2 求和与求积 数据序列求和与求积的函数是sum和prod，其使用方法类似。设X是一个向量，A是一个矩阵，函数的调用格式为： sum(X)：返回向量X各元素的和。 prod(X)：返回向量X各元素的乘积。 sum(A)：返回一个行向量，其第i个元素是A的第i列的元素和。 prod(A)：返回一个行向量，其第i个元素是A的第i列的元素乘积。 sum(A,dim)：当dim为1时，该函数等同于sum(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素之和。 prod(A,dim)：当dim为1时，该函数等同于prod(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素乘积。 例6-4 求矩阵A的每行元素的乘积和全部元素的乘积。 6.1.3 平均值和中值 求数据序列平均值的函数是mean，求数据序列中值的函数是median。两个函数的调用格式为： mean(X)：返回向量X的算术平均值。 median(X)：返回向量X的中值。

数值分析MATLAB上机实验

数值分析实习报告 姓名：gestepoA 学号：201******* 班级：***班

序言 随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。而且还能减少大量的人工计算。 由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JAVA的辅助软件来解决，得到一个满足误差限的解。本文所计算题目，均采用MATLAB进行编程，MATLAB被称为第四代计算机语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。它具有以下优点： 1友好的工作平台和编程环境。MATLAB界面精致，人机交互性强，操作简单。 2简单易用的程序语言。MATLAB是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（M 文件）后再一起运行。 3强大的科学计算机数据处理能力。包含大量计算算法的集合，拥有600多个工程中要用到的数学运算函数。 4出色的图像处理功能，可以方便地输出二维图像，便于我们绘制函数图像。

目录 1 第一题 (4) 1.1 实验目的 (4) 1.2 实验原理和方法 (4) 1.3 实验结果 (5) 1.3.1 最佳平方逼近法 (5) 1.3.2 拉格朗日插值法 (7) 1.3.3 对比 (8) 2 第二题 (9) 2.1实验目的 (9) 2.2 实验原理和方法 (10) 2.3 实验结果 (10) 2.3.1 第一问 (10) 2.3.2 第二问 (11) 2.3.3 第三问 (11) 3 第三题 (12) 3.1实验目的 (12) 3.2 实验原理和方法 (12) 3.3 实验结果 (12) 4 MATLAB程序 (14)

《MATLAB与数值分析》第一次上机实验报告

电子科技大学电子工程学院标准实验报告（实验）课程名称MATLAB与数值分析 学生姓名：李培睿 学号：2013020904026 指导教师：程建

一、实验名称 《MATLAB与数值分析》第一次上机实验 二、实验目的 1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算 操作。（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序） 2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号 转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。（用.m文件编写进行符号因式分解和函数求反的程序） 3. 掌握Matlab函数的编写规范。 4、掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、 三维曲线和面的填充、三维等高线等。（用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释） 5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。 三、实验内容 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以x， y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图，允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。不允许使用sort 函数。 四、实验数据及结果分析 题目一： ①在Editor窗口编写函数代码如下：

MATLAB与数值分析课程总结

MATLAB与数值分析课程总结 姓名：董建伟 学号：2015020904027 一：MATLAB部分 1.处理矩阵－容易 矩阵的创建 （1）直接创建注意 a中括号里可以用空格或者逗号将矩阵元素分开 b矩阵元素可以是任何MATLAB表达式，如实数复数等 c可以调用赋值过的任何变量，变量名不要重复，否则会被覆盖 （2）用MATLAB函数创建矩阵如：a空阵［］ b rand/randn——随机矩阵 c eye——单位矩阵 d zeros ——0矩阵 e ones——1矩阵 f magic——产生n阶幻方矩阵等 向量的生成 （1）用冒号生成向量 （2）使用linspace和logspace分别生成线性等分向量和对 数等分向量 矩阵的标识和引用 （1）向量标识 （2）“0 1”逻辑向量或矩阵标识 （3）全下标，单下标，逻辑矩阵方式引用 字符串数组 （1）字符串按行向量进行储存 （2）所有字符串用单引号括起来 （3）直接进行创建 矩阵运算 （1）注意与数组点乘，除与直接乘除的区别，数组为乘方对应元素的幂

（2）左右除时斜杠底部靠近谁谁是分母 （3）其他运算如，inv矩阵求逆，det行列式的值, eig特征值，diag 对角矩阵 2.绘图－轻松 plot－绘制二维曲线 （1）plot(x)绘制以x为纵坐标的二维曲线 plot(x,y) 绘制以x为横坐标，y为纵坐标的二维曲线 x,y为向量或矩阵 （2）plot（x1,y1,x2,y2，。。。。。。）绘制多条曲线，不同字母代替不同颜色:b蓝色，y黄色，r红色，g绿色 （3）hold on后面的pl ot图像叠加在一起 hold off解除hold on命令，plot将先冲去窗口已有图形（4）在hold后面加上figure，可以绘制多幅图形 （5）subplot在同一窗口画多个子图 三维图形的绘制 (1)plot3(x,y,z,’s’) s是指定线型，色彩，数据点形的字 符串 (2)［X,Y］＝meshgrid(x,y)生成平面网格点 (3)mesh(x,y,z，c)生成三维网格点，c为颜色矩阵 (4)三维表面处理mesh命令对网格着色，surf对网格片着色 (5)contour绘制二维等高线 (6)axis(［x1,xu,y1,yu］)定义x,y的显示范围 3．编程－简洁 （1）变量命名时可以由字母，数字，下划线，但是不得包含空格和标点 （2）最常用的数据类型只有双精度型和字符型，其他数据类型只在特殊条件下使用 （3）为得到高效代码，尽量提高代码的向量化程度，避免使用循环结构

matlab实现数值分析插值及积分

Matlab实现数值分析插值及积分 摘要： 数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。 分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。 其中Aitken插值计算的结果图如下： 对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932 其中复化梯形公式计算的结果图如下：

问题重述 问题一：已知列表函数 表格 1 分别用拉格朗日，牛顿，埃特金插值方法计算。 问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。 问题解决 问题一：插值方法 对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。 一、拉格朗日插值法： 拉格朗日插值多项式如下： 首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数 )(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式，所以只 可能相差一个常数因子，固)(x l i 可表示成： )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得：

数值分析的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程序

列主元法 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')

第6章matlab数据分析与多项式计算_习题答案

第6章 MATLAB数据分析与多项式计算 习题6 一、选择题 1．设A=[1,2,3,4,5;3,4,5,6,7]，则min(max(A))的值是（）。B A．1 B．3 C．5 D．7 2．已知a为3×3矩阵，则运行mean(a)命令是（）。B A．计算a每行的平均值 B．计算a每列的平均值 C．a增加一行平均值 D．a增加一列平均值 3．在MATLAB命令行窗口输入下列命令： >> x=[1,2,3,4]; >> y=polyval(x,1); 则y的值为（）。 D A．5 B．8 C．24 D．10 4．设P是多项式系数向量，A为方阵，则函数polyval(P,A)与函数polyvalm(P,A)的值（）。D A．一个是标量，一个是方阵 B．都是标量 C．值相等 D．值不相等 5．在MATLAB命令行窗口输入下列命令： >> A=[1,0,-2]; >> x=roots(A); 则x(1)的值为（）。 C A．1 B．-2 C． D． 6．关于数据插值与曲线拟合，下列说法不正确的是（）。A A．3次样条方法的插值结果肯定比线性插值方法精度高。 B．插值函数是必须满足原始数据点坐标，而拟合函数则是整体最接近原始数据点，而不一定要必须经过原始数据点。 C．曲线拟合常常采用最小二乘原理，即要求拟合函数与原始数据的均方误差达到极小。 D．插值和拟合都是通过已知数据集来求取未知点的函数值。 二、填空题 1．设A=[1,2,3;10 20 30;4 5 6]，则sum(A)= ，median(A)= 。 [15 27 39]，[4 5 6[ 2．向量[2,0,-1]所代表的多项式是。2x2-1 3．为了求ax2+bx+c=0的根，相应的命令是（假定a、b、c已经赋值）。为了

同济大学数值分析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实现 一、实验目的 学会在MATLAB环境下对已知的数据进行处理。 二、实验方法 1. 求取数据的最大值或最小值。 2. 求取向量的均值、标准方差和中间值。 3．在MATLAB环境下，对已知的数据分别进行曲线拟合和插值。 三、实验设备 1．586以上微机，16M以上内存，400M硬盘空间，2X CD-ROM 2．MATLAB5.3以上含CONTROL SYSTEM TOOLBOX。 四、实验内容 1．在MATLAB环境下，利用MATLAB控制系统工具箱中的函数直接求取数据的最大值或最小值，以及向量的均值、标准方差和中间值。 2．在MATLAB环境下，选择合适的曲线拟合和插值方法，编写程序，对已知的数据分别进行曲线拟合和插值。 五、实验步骤 1. 在MATLAB环境下，将已知的数据存到数据文件mydat.mat中。 双击打开Matlab，在命令窗口（command window）中，输入一组数据：实验一数据处理方法的MATLAB实现 一、实验目的 学会在MATLAB环境下对已知的数据进行处理。 二、实验方法 1. 求取数据的最大值或最小值。 2. 求取向量的均值、标准方差和中间值。 3．在MATLAB环境下，对已知的数据分别进行曲线拟合和插值。 三、实验设备 1．586以上微机，16M以上内存，400M硬盘空间，2X CD-ROM 2．MATLAB5.3以上含CONTROL SYSTEM TOOLBOX。 四、实验内容

1．在MATLAB环境下，利用MATLAB控制系统工具箱中的函数直接求取数据的最大值或最小值，以及向量的均值、标准方差和中间值。 2．在MATLAB环境下，选择合适的曲线拟合和插值方法，编写程序，对已知的数据分别进行曲线拟合和插值。 五、实验步骤 1. 在MATLAB环境下，将已知的数据存到数据文件mydat.mat中。 双击打开Matlab，在命令窗口（command window）中，输入一组数据： x=[1,4,2,81,23,45] x = 1 4 2 81 2 3 45 单击保存按钮，保存在Matlab指定目录（C:\Program Files\MATLAB71）下，文件名为“mydat.mat”。 2. 在MATLAB环境下，利用MATLAB控制系统工具箱中的函数直接求取数据的最大值或最小值，以及向量的均值、标准方差和中间值。 继续在命令窗口中输入命令： （1）求取最大值“max（a）”； >> max(x) ans = 81 （2）求取最小值“min（a）”； >> min(x) ans = 1 （3）求取均值“mean（a）”； >> mean(x) ans =

数值分析算法在matlab中的实现

数值分析matlab实现高斯消元法： function[RA,RB,n,X]=gaus(A,b) B=[A b];n=length(b);RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1 for k=p+1:n m=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else disp('请注意：因为RA=RB0, disp('请注意：因为RA~=RB，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1

数值分析matlab函数资料

1．求数值积分： fx=@(x)exp(1./x); >> quadl(fx,1,5) 2．获取x=xlsread('oillack.xls','sheet1','a1:a73') excel文件名是oillack.xls，sheet1是表名，a1:a73'是a列的1到73行 long x=xlsread('F:\学习\大三\大三下\巷道力学模型\新建文件夹(2)\1.xlsx','sheet1','a2:a') 3. 在matlab的图中插入文本框后将文本框旋转的方法： text(0.5,0.6,'渗透率/mD','Rotation',90) 4. matlab中插入一条直线的方法： line([0.01 0.01],[0 1.75]) 5．Matlab 中画三维图 x=-7.5:0.5:7.5; y=x; % 先产生x及y二个阵列 >> [x,y]=meshgrid(x,y); % 再以meshgrid形成二维的网格数据 >> z=x.^2+y.^2; % 产生z轴的数据 >> mesh(x,y,z) % 将z轴的变化值以网格方式画出 >> surf(x,y,z) % 将z轴的变化值以曲面方式画出 Matlab指数拟合方法 x=[1982 1992 2002]； y=[103.5 34.5 23.3]； cftool（x，y） 在弹出的对话框选择fitting，弹出新的对话框选择new fit，然后在第三个下拉菜单（Type of fit）中选择Exponential，然后点击Apply，即可；最后结果 General model Exp1: f(x) = a*exp(b*x) Coefficients (with 95% confidence bounds): a = 1.453e+082 (-7.288e+084, 7.317e+084) b = -0.09312 (-0.3464, 0.1602)

数值分析matlab代码

1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6) disp('牛顿法'); i=1; n0=180; p0=pi/3; tol=10^(-6); for i=1:n0 p=p0-(p0-sin(p0))/(1-cos(p0)); if abs(p-p0)<=10^(-6) disp('用牛顿法求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次牛顿迭代后无法求出方程的解') end 2、disp('Steffensen加速'); p0=pi/3; for i=1:n0 p1=0.5*p0+0.5*cos(p0); p2=0.5*p1+0.5*cos(p1); p=p0-((p1-p0).^2)./(p2-2.*p1+p0); if abs(p-p0)<=10^(-6) disp('用Steffensen加速求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次Steffensen加速后无法求出方程的解') end 1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根， %误差限为 e=10^-4 disp('二分法')

a=0.2;b=0.26; tol=0.0001; n0=10; fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1; for i=1:n0 p=(a+b)/2; fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1; if fp==0||(abs((b-a)/2)0 a=p; else b=p; end end if i==n0&&~(fp==0||(abs((b-a)/2)

数值分析 matlab 实验4

(1) 解题过程如下： （1）MATLAB中创建复化梯形公式和复化辛普森公式的 M 文件：1）复化梯形公式文件： function s=T_fuhua(f,a,b,n) h=(b-a)/n; s=0; for k=1:(n-1) x=a+h*k; s=s+feval(f,x); end s=h*(feval(f,a)+feval(f,b))/2+h*s; 2）复化辛普森公式文件： function s=S_fuhua(f,a,b,n) h=0; h=(b-a)./(2*n); s1=0; https://www.wendangku.net/doc/9814967178.html, -5- s2=0; for k=1:n-1 x=a+h*2*k; s1=s1+feval(f,x); end for k=1:n x=a+h*(2*k-1); s2=s2+feval(f,x); end

s=h*(feval(f,a)+feval(f,b)+s1*2+s2*4)/3; 在MATLAB中输入： f=inline('x/(4+x^2)');a=0;b=1; %inline 构造内联函数对象 for n=2:10 s(n-1)=T_fuhua(f,a,b,n);s(n-1)=vpa(s(n-1),10); %调用复化梯形公式，生成任意精度的数值 end exact=int('x/(4+x^2)',0,1);exact=vpa(exact,10) %求出积分的精确值 输出结果：exact = .1115717755 s = Columns 1 through 6 0.1088 0.1104 0.1109 0.1111 0.1113 0.1114 Columns 7 through 9 0.1114 0.1114 0.1115 在MATLAB中输入以下函数用以画出计算误差与 n 之间的曲线： r=abs(exact-s); n=2:10; plot(double(n),double(r(n-1))) 得到结果如图所示： （2）在 MATLAB中输入以下程序代码： f=inline('x/(4+x^2)');a=0;b=1;n=9; %inline 构造内联函数对象 t=T_fuhua(f,a,b,n);t=vpa(t,10) s=S_fuhua(f,a,b,n);s=vpa(s,10)

数值分析实验— MATLAB实现

数值分析实验 ——MATLAB实现 姓名sumnat 学号2013326600000 班级13级应用数学2班 指导老师 2016年1月

一、插值：拉格朗日插值 (1) 1、代码： (1) 2、示例： (1) 二、函数逼近：最佳平方逼近 (2) 1、代码： (2) 2、示例： (2) 三、数值积分：非反常积分的Romberg算法 (3) 1、代码： (3) 2、示例： (4) 四、数值微分：5点法 (5) 1、代码： (5) 2、示例： (6) 五、常微分方程：四阶龙格库塔及Adams加速法 (6) 1、代码：四阶龙格库塔 (6) 2、示例： (7) 3、代码：Adams加速法 (7) 4、示例： (8) 六、方程求根：Aitken 迭代 (8) 1、代码： (8) 2、示例： (9) 七、线性方程组直接法：三角分解 (9) 1、代码： (9) 2、示例： (10) 八、线性方程组迭代法：Jacobi法及G-S法 (11) 1、代码：Jacobi法 (11) 2、示例： (12) 3、代码：G-S法 (12) 4、示例： (12) 九、矩阵的特征值及特征向量：幂法 (13) 1、代码： (13) 2、示例： (13)

一、插值：拉格朗日插值 1、代码： function z=LGIP(x,y)%拉格朗日插值 n=size(x); n=n(2);%计算点的个数 syms a; u=0;%拉格朗日多项式 f=1;%插值基函数 for i=1:n for j=1:n if j==i f=f; else f=f*(a-x(j))/(x(i)-x(j)); end end u=u+y(i)*f;f=1; end z=expand(u);%展开 2、示例： >> x=1:6; y1=x.^5+3*x.^2-6; y2=sin(x)+sqrt(x); >> f1=LGIP(x,y1) f1 = -6+3*a^2+a^5 %可知多项式吻合得很好 >> f2=vpa(LGIP(x,y2),3) f2 = .962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5

数值分析重要算法的matlab程序

数值分析重要算法的matlab程序 插值多项式:拉格朗日插值 function yh=lagrange(x,y,xh) n=length(x); m=length(xh); yh=zeros(1,m); c1=ones(n-1,1); c2=ones(1,m); for i=1:n xp=x([1:i-1 i+1:n]); yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2)); end 插值多项式:牛顿插值（以数值实验3.2）为例 function y=ex32 n = 21; x = linspace(-5,5,n)'; h = (5-(-5))/(n-1); y = 1./(1+x.^2); % form the differences table for j = 2:n, y(1:n+1-j,j) = diff(y(1:n+2-j,j-1))./(x(j:n)-x(1:n+1-j)); end % newton coeff y = y(1,:); pz = [ ]; v = linspace(-5,5,80); for t = v, z = y(n); for j = n-1:-1:1, z = z * (t-x(j))+y(j); end pz=[pz z]; end plot(v,pz,'r+-',v,1./(1+v.^2),'g--'); 数值积分：梯形求积公式求积分 function I=ftrapz(fun,a,b,n) h=(b-a)/n; x=linspace(a,b,n+1); y=feval(fun,x); I=h*(0.5*y(1)+sum(y(2:n))+0.5*y(n+1)); 数值积分：抛物型求积公式求积分 function I=fsimpsion(fun,a,b,n) h=(b-a)/n; x=linspace(a,b,2*n+1); y=feval(fun,x); I=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1)); 追赶法解三对角方程组

数值计算方法matlab程序

function [x0,k]=bisect1(fun1,a,b,ep) if nargin<4 ep=1e-5; end fa=feval(fun1,a); fb=feval(fun1,b); if fa*fb>0 x0=[fa,fb]; k=0; return; end k=1; while abs(b-a)/2>ep x=(a+b)/2; fx=feval(fun1,x); if fx*fa<0 b=x; fb=fx; else a=x; fa=fx;

end end x0=(a+b)/2; >> fun1=inline('x^3-x-1'); >> [x0,k]=bisect1(fun1,1.3,1.4,1e-4) x0 = 1.3247 k = 7 >> 简单迭代法 function [x0,k]=iterate1(fun1,x0,ep,N) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep;

while abs(x-x0)>ep & k> fun1=inline('(x+1)^(1/3)'); >> [x0,k]=iterate1(fun1,1.5) x0 = 1.3247 k = 7 >> fun1=inline('x^3-1'); >> [x0,k]=iterate1(fun1,1.5) x0 = Inf k =

基于MATLAB《数值分析》有关算法的实现

基于MATLAB 《数值分析》有关算法的实现 摘要： 数值分析算法在科学应用各个分支中有着越来越广泛的应用，特别是计算机的出现更有助于这些应用的不断发展。实际生活中，很多问题都可以归纳为数值问题，如测得的实验数据有时候可以用插值法进行函数逼近预测等。数值算法常与工程实践相结合，通过MATLAB 软件运行的程序来解决实际问题。本课题实现了拉格朗日插值法、最小二乘法、雅克比迭代法、二分法等算法的MATLAB 程序设计。 拉格朗日插值法 1. 基本原理： 通过平面上不同的两点可以确定一条直线，这就是拉格朗日线性插值问题，对于不在同一条直线的三个点得到的插值多项式则为抛物线。 拉格朗日插值的基本多项式为： l i (x )=∏x?x i x i ?x j j=0j≠i ,i=0,1,2,…n 有了基函数以后就可以直接构造如下多项式： p n =∑f(x i )l i (x)n i=0 该多项式就是拉格朗日插值法所求得的插值多项式。 2. 拉格朗日插值法算法： 1、 根据所给点（x i ，y i ）的坐标依次写出其插值基函数（用for 循环可以轻易解决） l i (x )=∏x?x i x i ?y i ,j=0j≠i i=0,1,2,…,n 2、 将插值基函数与其对应的点的函数值相乘得： f(x i )l i (x),i=0,1,2,…,n 3、 将2中各项累加即得插值多项式： P n =∑f(x i )l i (x)n i=0 3. 拉格朗日插值法程序： function lagrange(A) %A 为一个只有两行的矩阵，第一行为插值点，第二行 为插值点对应的函数值

数值分析matlab程序

数值分析（matlab程序） 曹德欣曹璎珞 第一章绪论 数值稳定性程序，计算P11 试验题一积分 function try_stable global n global a a=input('a='); N = 20; I0 = log(1+a)-log(a); I = zeros(N,1); I(1) = -a*I0+1; for k = 2:N I(k) = - a*I(k-1)+1/k; end II = zeros(N,1); if a>=N/(N+1) II(N) = (1+2*a)/(2*a*(a+1)*(N+1)); else II(N) =(1/((a+1)*(N+1))+1/N)/2; end for k = N:-1:2 II(k-1) = ( - II(k) +1/k) / a; end III = zeros(N,1); for k = 1:N n = k; III(k) = quadl(@f,0,1); end clc fprintf('\n 算法1结果算法2结果精确值') for k = 1:N, fprintf('\nI(%2.0f) %17.7f %17.7f %17.7f',k,I(k),II(k),III(k)) end function y = f(x) global n global a y = x.^n./(a+x); return

第二章非线性方程求解 下面均以方程y=x^4+2*x^2-x-3为例： 1、二分法 function y=erfen(a,b,esp) format long if nargin<3 esp=1.0e-4; end if fun(a)*fun(b)<0 n=1; c=(a+b)/2; while c>esp if fun(a)*fun(c)<0 b=c; c=(a+b)/2; elseif fun(c)*fun(b)<0 a=c; c=(a+b)/2; else y=c; esp=10000; end n=n+1; end y=c; elseif fun(a)==0 y=a; elseif fun(b)==0 y=b; else disp('these,nay not be a root in the intercal') end n function y=fun(x) y=x^4+2*x^2-x-3; 2、牛顿法 function y=newton(x0) x1=x0-fun(x0)/dfun(x0); n=1; while (abs(x1-x0)>=1.0e-4) & (n<=100000000) x0=x1; x1=x0-fun(x0)/dfun(x0); n=n+1; end y=x1 n function y=fun(x)