紧急求助:三次样条逼近的程序

Python017

紧急求助:三次样条逼近的程序,第1张

我们学校的数学建模上机课也有Mathlab程序,看看下面有没有你要找的。

一 基本运算

1 求

输入(12+2*(7-4))/3^2执行

2 输入x = (5*2+1.3-0.8)*10^2/25执行

再输入y= 2*x+1执行

3 执行clear命令。观察结果

4计算圆面积Area = ,半径r = 2,则可键入

r=2area=pi*r^2area

问:语句末尾加分号与不加分号有何区别?请试验之

5常用函数

名称 含义 名称 含义

sin 正弦 exp E为底的指数

cos 余弦 log 自然对数

tan 正切 log10 10为底的对数

cot 余切 log2 2为底的对数

asin 反正弦 abs 绝对值

acos 反余弦

例:1)执行y = sin(10)*exp(-0.3*4^2)

2) 想计算 的值

输入y1=2*sin(0.3*pi)/(1+sqrt(5))执行之

若又想计算 ,可以简便地用操作:先按á键则会出现上面输入过的指令 y1=2*sin(0.3*pi)/(1+sqrt(5)) ;然后移动光标,把y1改成y2;把 sin 改成 cos 便可。即得

y2=2*cos(0.3*pi)/(1+sqrt(5))然后执行之。

系统默认4位有效数字,若想提高精度则可如下:

digits(10)sym(y2,'d') 执行就可精确到小数点后10位,还可将10改为其它数字试验

矩阵运算

1要得到矩阵 ,

可输入A = [1,2,34,5,67,8,9] 执行,观察结果

还可分行输入

A=[1,2,3

4,5,6

7,8,9]

效果相同

2 注意 %号后的语句为注释,练习时不必输入

>>a=[1,4,6,8,10] %一维矩阵

>>a(3) % a的第三个元素

ans =

6

»x =[1 2 3 4 5 6 7 8

4 5 6 7 8 9 10 11]%二维2x8 矩阵

执行后双击左边Workspace里的x,观察之

» x(3) % x的第三个元素

ans =

2

» x([1 2 5]) % x的第一、二、五个元素

ans =

1 4 3

如需要还可定义b=x([1 2 5])执行后结果为

b =

1 4 3

>>x(2,3) % x的第二行第三列的元素

ans =

6

x(1:5) % x的第前五个元素

ans =

1 4 2 5 3

» x(10:end) % x的第十个元素后的元素

ans =

8 6 9 7 10 8 11

执行后双击左边Workspace里的x,观察是哪十个元素

» x(10:-1:2) % x的第十个元素和第二个元素的倒排

ans =

8 5 7 4 6 3 5 2 4

» x(find(x>5)) % x中大于5的元素

ans =

6 7 8 6 9 7 10 8 11

» x(4)=100 %给x的第四个元素重新给值

x =

1 2 3 4 5 6 7 8

4 100 6 7 8 9 10 11

» x(3)=[] % 删除第三个元素(不是二维数组)

x =

Columns 1 through 12

1 4 100 3 6 4 7 5 8 6 9 7

Columns 13 through 15

10 8 11

» x(16)=1 % 加入第十六个元素

x =

Columns 1 through 12

1 4 100 3 6 4 7 5 8 6 9 7

Columns 13 through 16

10 8 11 1

3 如不需要以前的变量时,为不干扰以后计算,可执行clear清除以前的变量

当元素很多的时候,则须采用以下的方式:

» x=(1:2:121)% 以起始值为1,增量值为2,终止值为121的矩阵

» x=linspace(0,1,100)% 利用linspace,生成以0为起始值,1为终止值,元素数目为100的矩阵

»a=[] %空矩阵

a =

[]

» zeros(2,2) %全为0的矩阵

ans =

0 0

0 0

» ones(3,3) %全为1的矩阵

ans =

1 1 1

1 1 1

1 1 1

» rand(2,4)%随机矩阵

4另外一种定义矩阵的方式

»a=1:7b=11:2:23

»c=[b a]%利用上面建立的阵列 a 及阵列 b ,组成新阵列c

»d=[b a]%利用a及b,组成新矩阵d

执行后双击左边Workspace里的c与d,比较之

再如 已知y=[-1,6,15,7,31,2,4,5]

x=y(3:5) %x为y的第三到第五个元素组成的新向量

或 x=[y(5),y(3),y(7)] %x为y的第五、第三、第七个元素组成的新向量

或这样更简单 x=y([5,3,7])

5 输入矩阵x=[4,8,12,10,236,3,15,13,199,1,2,18,1411,7,5,21,17]

依次输入下列命令并执行,观察结果,各命令分别有什么作用?

max(x)

min(x) (问:如何得到整个矩阵的最小值与最大值?)

[m,n]=size(x)

L=length(x)

y=x’

a=x( :,2)

b=x( :,2)’

c=x(3, :)

d=x(1 :3,3 :5)

y(2,3)=y(2,3)/2

y(2, :)=y(2, :)/2

y( :,4)=y( :,4)+y( :,2)

6 点运算 执行下列命令,指出点运算的作用

x=1 :8 (或对另外的向量或矩阵来作)

y=2.^x

z=x./y

w=x.^2

u=sin(x)

常用命令

min 最小值 max 最大值

mean 平均值 std 标准差

sort 排序 diff 相邻元素的差

length 个数 sum 总和

dot 内积 cross 外积

三 画图

二维图形

命 令 含 义 plot绘图函数的叁数

plot 建立向量或矩阵各队队向量的图形 字元 颜色 字元 图线型态

loglog x、y轴都取对数标度建立图形 y 黄色 . 点

semilogx x轴用于对数标度,y轴线性标度绘制图形 k 黑色 o 圆

semilogy y轴用于对数标度,x轴线性标度绘制图形 w 白色 x x

title 给图形加标题 b 蓝色 + +

xlabel 给x轴加标记 g 绿色 * *

ylabel 给y轴加标记 r 红色 - 实线

text 在图形指定的位置上加文本字符串 c 亮青色 : 点线

gtext 在鼠标的位置上加文本字符串 m 锰紫色 -. 点虚线

grid 打开网格线 -- 虚线

hold on 命令用于在已画好的图形上添加新的图形

1 x=0:0.001:10% 0到10的1000个点(每隔0.001画一个点)的x座标

y=sin(x)% 对应的y座标

plot(x,y)% 绘图

注:matlab画图实际上就是描点连线,因此如果点取得不密,画出来就成了折线图,请试验之

2 Y=sin(10*x)

plot(x,y,'r:',x,Y,'b') % 同时画两个函数

3 若要改变颜色,在座标对后面加上相关字串即可:

x=0:0.01:10

plot(x,sin(x),'r')

4 若要同时改变颜色及图线型态(Line style),也是在坐标对后面加上相关字串即可:

plot(x,sin(x),'r*')

5 用axis([xmin,xmax,ymin,ymax])函数来调整图轴的范围

axis([0,6,-1.5,1])

6 MATLAB也可对图形加上各种注解与处理:(见上表)

xlabel('x轴')% x轴注解

ylabel('y轴')% y轴注解

title('余弦函数')% 图形标题

legend('y = cos(x)')% 图形注解

gtext('y = cos(x)')% 图形注解 ,用鼠标定位注解位置

grid on% 显示格线

7画椭圆

a = [0:pi/50:2*pi]'%角度

X = cos(a)*3%参数方程

Y = sin(a)*2

plot(X,Y)

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

title('椭圆')

8 绘制函数 在0 ≤ x ≤ 1时的曲线。

x=0:0.1:1

y=x.*exp(-x) %为什么用点运算?若不用会怎样

plot(x,y),xlabel('x'),ylabel('y'),title('y=x*exp(-x)')

9 画出衰减振荡曲线 与它的包络线 及 。t 的取值范围是[0, 4π] 。

t=0:pi/50:4*pi

y0=exp(-t/3)

y=exp(-t/3).*sin(3*t)

plot(t,y,'-r',t,y0,':b',t,-y0,':b') % -r表示红色实线,:b表示蓝色点线,看上表

grid

10 在同一个画面上建立几个坐标系, 用subplot(m,n,p)命令;把一个画面分成m×n个图形区域, p代表当前的区域号,在每个区域中分别画一个图,如

x=linspace(0,2*pi,30)y=sin(x)z=cos(x)

u=2*sin(x).*cos(x)v=sin(x)./cos(x)

subplot(2,2,1),plot(x,y),axis([0 2*pi -1 1]),title('sin(x)')

subplot(2,2,2),plot(x,z),axis([0 2*pi -1 1]),title('cos(x)')

subplot(2,2,3),plot(x,u),axis([0 2*pi -1 1]),title('2sin(x)cos(x)')

subplot(2,2,4),plot(x,v),axis([0 2*pi -20 20]),title('sin(x)/cos(x)')

三维图形

11三维螺旋线:

t=0:pi/50:10*pi

plot3(sin(t),cos(t),t) %参数方程

grid %添加网格

12 t=linspace(0,20*pi, 501)

plot3(t.*sin(t), t.*cos(t), t)%注意点乘

也可以同时画出两条曲线,格式与二维情况类似,兹不举例。

13用mesh命令画曲面

画出由函数 形成的立体网状图:

a=linspace(-2, 2, 25)% 在x轴上从(-2,2)取25点

b=linspace(-2, 2, 25)% 在y轴上取25点

[x,y]=meshgrid(a, b)% x和y都是21x21的矩阵

z=x.*exp(-x.^2-y.^2)% 计算函数值,z也是21x21的矩阵

mesh(x, y, z)% 画出立体网状图

14 surf和mesh的用法类似:

a=linspace(-2, 2, 25)% 在x轴上取25点

b=linspace(-2, 2, 25)% 在y轴上取25点

[x,y]=meshgrid(a, b)% x和y都是21x21的矩阵

z=x.*exp(-x.^2-y.^2)% 计算函数值,z也是21x21的矩阵

surf(x, y, z)% 画出立体曲面图

四 程序设计

1 M-文件: 上面所做的运算都是在命令窗口中输入一条或两三条命令,然后执行,再输入,再执行,以这样交谈式的方式进行。如果为了解决某一问题需要很多命令,这样做就很不方便了。这时我们把解决某一问题的所有命令集中放在一个文档里,命名、保存。然后只要在命令窗口中输入文档名,执行即可。

例:(1)编写文档:点击MATLAB指令窗口上面最左端的图标 ,即新建文件,就可打开MATLAB文件编辑器。用户即可在空白窗口中编写程序。例如输入下面的程序:

x=linspace(0,2*pi,20)

y=sin(x)

plot(x,y,'r+')

title('2D plot')

(2)点击文件编辑器上面工具条中的保存 ,命名(例如将上面的程序命名为picture),然后保存。像这样在MATLAB文件编辑器中编写的文件叫M-文件(M-file)。

(3)运行:i)在命令窗口中输入文件名(如上面的picture),然后执行。

ii)或直接在文件编辑器上面的工具条中找到debug(即调试),点击,再找到run(即运行),再点击即可。

同学们可以把前面画图的一些问题放在文件编辑器里再做一下。

2 自己编写函数:我们经常用到的像sin、cos、exp这样的一些函数都是MATLAB软件自身所带的函数,因此直接应用即可,但有时我们为了解决一些问题需要自己编写函数。自己编写函数有两个基本要求i)必须在MATLAB文件编辑器中编写。ii)函数名和文件名必须相同。 例: 编写函数 , 计算f(1)f(2)+f2(3)

(1)打开MATLAB文件编辑器,输入

function Y= fun1(x) % 表示Y是x的函数,x是自变量,fun1是函数名

Y=(x^3 - 2*x^2 + x - 6.3)/(x^2 + 0.05*x - 3.14)

然后保存。

注:在自己编写的函数前都要写上function,表示这是自己定义的函数。fun1表示函数名,那么最后文件名也应命名为fun1。

(2)这样在命令窗口中就可以像应用sin、cos那样来使用函数fun1,如:在命令窗口中输入 >>fun1(1)*fun1(2)+fun1(3)*fun1(3) 结果为:

ans =

-12.6023

3 for循环语句(这里的for语句与C语言中的for语句不同,要更简单一些)

例:一个简单的for循环示例。

for i=1:10% i依次取1,2,…10,.

x(i)=2*i% 对每个i值,重复执行该指令

end% 表示循环结束,每一个for要对应一个end

x % 要求显示运行后数组x的值。

输入后观察结果,体会for语句的作用。

注:在MATLAB里(在C语言中也一样), 的作用表示把等号右边的值送给左边的变量,这和数学中相等的意思不同。下面的例子中都要这样理解,否则就不能明白程序的含义。

4 while循环语句

例: Fibonacci 数列:1,1,2,3,5,8,… 即: ,( 1,2,3…)现要求该数列中第一个大于10000 的元素。

a(1)=1a(2)=1i=2

while a(i)<=10000

a(i+1)=a(i-1)+a(i)

i=i+1

end

i,a(i),

5(1)if-end语句,例:

cost=10number=12

if number>8

sums=number*0.95*cost

end,

sums

(2)if-else-end语句,例:

cost=10number=5% 改变number的初值,看结果有何不同

if number>8

sums=number*0.95*cost

else sums=number*0.5*cost

end,

sums

6 例:用for 循环语句来寻找Fibonacc 数列中第一个大于10000 的元素。

n=100a=ones(1,n)% a是一个一行,n列的所有元素为1的矩阵

for i=3:n

a(i)=a(i-1)+a(i-2)

if a(i)>=10000

a(i),

break% 表示跳出循环

end

end, i

7 练习:课本264页,参考例4右边的流程图11.4,编程序求解例4,自己设置误差,并与书上的结果比较。

五 拟合与插值

曲 线 拟 合 和 插 值 函 数

polyfit(x, y, n) 对描述n阶多项式y=f(x)的数据进行最小二乘曲线拟合

interp1(x, y, xo) 1维线性插值

interp1(x, y, xo, ' spline ') 1维3次样条插值

interp1(x, y, xo, ' cubic ') 1维3次插值

interp2(x, y, Z, xi, yi) 2维线性插值

interp2(x, y, Z, xi, yi, ' cubic') 2维3次插值

1 插值

看课本266页§11.2第一段,了解什么是插值。

例:考虑下列问题,12小时内,一小时测量一次室外温度。数据如下:

时间:1,2,3, 4, 5, 6, 7, 8, 9,10,11,12

温度:5,8,9,15,25,29,31,30,22,25,27,24

现在根据以上数据估计3.2,4.7等时刻的温度

hours=1:12

temps=[5 8 9 15 25 29 31 30 22 25 27 24]

t=interp1(hours, temps, [3.2,4.7]) % 一阶线性插值,如果只估计一个点的值,则无须加方括号

改为t=interp1(hours, temps, [3.2,4.7], 'spline') 则为三次样条插值

如果输入如下程序,则画出插值曲线

hours=1:12

temps=[5 8 9 15 25 29 31 30 22 25 27 24]

h=1:0.1:12

t=interp1(hours, temps, h) % h后加上'spline'则为三次样条插值

plot(hours, temps, ' + ' , h, t)

用一阶线性插值和三次样条插值做课本268页例2,与书上之结果比较,然后挑课后题做一两道。

2 拟合

看课本270页§11.3,曲线拟合,比较拟合与插值有什么区别。

例:两组数据如下:

x=[0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1]

y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]

n=8

p1=polyfit(x,y,n)% n表示用n阶多项式拟合,n=1为线性拟合,即通常所说最小二乘法

poly2sym(p1) % 前面的拟合命令只给出多项式的系数,用此命令则将结果转化为真正的多项式。或用 vpa(poly2sym(p1),10) 即取数值形式,取10位有效数字

x1=0:.01:1% 由此以后三句是画出拟合曲线的图像

y1=polyval(p1,x1)%此句是在x1这些点处求出多项式的值,送给y1

plot(x,y,'o',x1,y1)

改变n的数字,即用不同的多项式拟合,看看哪个结果好。

当n=10时,数据点之间出现大的波动。当企图进行高阶曲线拟合时,这种波动现象经常发生,并不利于我们认识两组数据之间的规律,因此并不是阶数越高越好,实际问题当中,适当选一个即可。

用上面的指令做课本271页例1及例2,将结果与书上之结果比较。

2=6.610

x3=6.082

x4=5.607

x5=5.423

xH=5.825

%我这里的版本是2006a的,不知道6.5是否可以通过,你可以自己试试。

x=[x1 x2 x3 x4 x5]

y=[638.3 640.2 650.6 653.2 659.9]%你原来x只有5个元素,而y却有7个元素,所以会出错,我删了最后两个

xi=linspace(x1,x5,100)%这里自动让它平均生成了100个点,如果觉得太多可以自行改少一点

%你下面左引号都是全角中文的,应该换成英文半角状态

ycubic=interp1(x,y,xi,'cubic')

ylinear=interp1(x,y,xi,'linear')

yH1=interp1(x,y,xH,'cubic')

yH2=interp1(x,y,xi,'linear')

plot(x,y,'*',xH,yH1,'o',xi,ycubic,'-')

grid on

————————————————————————————

你们老师果然不负责任啊,选课的时候也不说明的。

建议自己找本matlab的教材学习一下,应该还是很好上手的。