用scilab通过调用脚本计算给定u和i的B样条基函数
由给定的u值来求B样条曲线上面的点的公式是:,Pi为控制顶点,U={u0,u1,……,un,}(这里的ui<uj,i<j),若u恰好在ui之后,C(u)只在控制顶点P(i-p),...,Pi处对应的基函数的值非零,所以上面的式子简化后改写的代码为
_3D_point Cu={0,0,0};
for(j=i-p;j<i;j++)
{
Cu=Cu+Nj,p(u)*Pj;
}
下面就给出求解B样条基的scilab脚本表示:
//文件名:nij.sce
//**********定义基函数**********通过给定的u求其函数值
function y=B_spline(i,u,p,U,N)
N(1)=1;
left=zeros(1,1);
right=zeros(1,1);
for j=1:p
left(j+1)=u-U(i+1-j);
right(j+1)=U(i+j)-u;
saved=0;
for r=0:(j-1)
temp=N(r+1)/(right(r+2)+left(j-r+1));
N(r+1)=saved+right(r+2)*temp;
saved=left(j-r+1)*temp;
end;
N(j+1)=saved;
disp(N); //显示每一次的N的值,在后面的输出结果中就是这里的显示end;
y=N;
endfunction
下面,我在scilab中对该函数进行调用,这里的U={0,0,0,1,2,3,4,4,5,5,5},控制顶点为随机的11个点,u=2.5,i=4,p=3
在scilab中输入下面的脚本,进行调用:
-->i=5;
-->u=2.5;
-->p=3;
-->U=[0,0,0,1,2,3,4,4,5,5,5];
-->exec('nij.sce');
-->y=zeros(1,1);
-->N=zeros(1,1);
-->y=B_spline(i,u,p,U,N); //下面是显示的结果(不包括j=0时 的结果):
0.5
0.5
0.125
0.75
0.125
0.0208333
0.4791667
0.46875
0.03125
1. %i^2 //-1
2. %nan==>not-a-namber
3. 如果你想输入一个3×4的向量,你就可以:a=[1 4 7 8;5 8 9 6;7 5 6 2],很简单。当然,你也可以这样:
a=[1 4 7 8
5 8 9 6
7 5 6 2]向量的转置:用单引号(‘),如a’就是将a转置。
4. s=poly(0,”cosQ”)或者s=poly(0,’cosQ’),s=poly(0,'s'); //poly內定为使用根來描述多項式 显示为s=cosQ
s^(1:10) 输出cosQ的从1到10 次方
s^1:10 输出cosQ:1:10
5. rand()用于产生随机数,如rand(2,3)产生随机的2 ×3矩阵
0.2113249 0.0002211 0.6653811
0.7560439 0.3303271 0.6283918
但是不知道为什么,数都是0~1之间的数
6. det(A)用于求矩阵的行列式
7. x=Ab 用于解线性方程组 A*x=b
8. [l u]=lu(A) 用于求矩阵A的lu分解,但是不知道为什么,分解的感觉不正确;
9. [l u e]=lu(A) 不太懂的lu分解,其中A(m×n),l(m×min(m,n)), u(min(m,n) ×n),e(n×n) l是下三角矩阵,u是上三角矩阵,e*A=l*u
10. l(i,j) 用于求矩阵l中的第i行,第j列元素,注意,这里的i、j均是从1开始的,而不是从0开始的,和c++语言中不一样;
11. sp=sparse([1,2;5,4;3,1],[1,2,3]) 用于构造稀疏矩阵,第1行第2列元素为:1,第5行第4列元素为:2,第3行第1列元素为:3
12. 查看全局变量用who(‘global’),清除全局变量b用clearglobal b
13. j:k表示[j,j+1,……,k];而j:d:k表示从d开始,每步增加d,直到小于等于k时为止;如一些具有固定步长量的向量如:[0 0.5 1 1.5 2 2.5 3] 可以用0:0.5:3来表示
14. A(:)表示显示矩阵A的所有元素,按列的顺序进行显示;
A(:,j)表示第j列元素
A(j:k)表示 [A(j),A(j+1),...,A(k)]
A(:,j:k)is [A(:,j),A(:,j+1),...,A(:,k)]
15. 求A的逆矩阵,在scilab中用inv(A)表示,求行列式的值,det(B),当然这里的B必须是个n*n矩阵。
16. 向量的点乘用 .*,如c=a.*b,./也一样,如a=[1,5,6],b=[2,2,5],c=a.*b结果为c=[2,10,30],a=c./b,结果为a=[1,5,6]
17. sum(A), 一个矩阵中所有数加在一起等于多少,用sum(X)来运算,
18. prod(A),矩阵的各项相乘
19. diag(A),提取对角元素
20. sqrt(平方根), cos(余弦),sin(正弦),max(最大值),round(四舍五入),括号里面都可以是数或者矩阵;
21. sign(判断正负),负值返回-1,正值返回1;
22. fft(a,1)快速傅立叶变换,
23. sort(a),递减排序,矩阵以列为主序号;
24. [a1,k]=gsort(a,‘g’,‘i’);
25. A=find(a>5),返回的是先列后行的顺序的元素序号,而不是具体的值;
26. zeros(m,n),生成m*n的零矩阵;eye(m*n),生成m*n的对角线元素为1的矩阵;ones(m*n)生成m*n 的由1组成的矩阵;
27. matrix(重组矩阵),感觉像是转值
28. det(求行列式值),inv(矩阵的求逆),qr(二次余数分解),svd(奇异值分解),bdiag(求广义本征值),spec(求本征值),schur(schur分解),trace(求对角线元素总和)
29. poly(构造多项式),如p=poly([1 2 3],'x')可能就是以1,2,3为根构造多项式,p=poly([-9 6 -5 4 1],'z','coeff') //p是z的多項式系数为,[ ] 内 0次是1次,……的系数,p= - 9 + 6z - 5z^2 + 4z^3 + z^4;roots(多项式求根),如roots(p)。由根来定义多项式方程式,poly([-1 -2+2*%i -2-2*%i -5+7*%i -5-7*%i],'z','root')。horner(y,2) //x=2時,y的数值。M=[p p-1; p+1 2] //多項式矩阵,det(M) //多项式矩阵M的行列式值。
30. coeff(求矩阵多项式的系数),horner(多项式求值),freq(求响应频率),clean
31. 文件加载运算,exec 用法-->exec("****.sce"),调试:pause,return,abort。样条函数,内插法:splin,interp,interpln。字符串:string,part,evstr,execstr
32. 绘图函数:plot,xset,driver,plot2d,xgrid,locate,plot3d,Graphics,要查询它们的用法,只需:help plot 就可以了,一般的说明文档都说的很清楚,而且还带有例子,只要感兴趣时去实践一下,掌握不成问题。
33. function y=fct1(x),y=x+exp(x)-10,endfunction //定义函数fct1,fsolve(1,fct1) //求解fct1=0的解,1为起始值
34. 微分方程式,dy/dt=y^2-y sin(t)+cos(t), y(0)=0, function ydot=f(t,y), ydot=y^2-y*sin(t)+cos(t),endfunction //定义微分方程式dy/dt=f(t,y);t0=0;y0=0; //起始值条件; t=0:0.1:%pi; //计算范围t由0开始到%pi(=3.1415927)间隔为0.1, y=ode(y0,t0,t,f); //呼叫ode函数求解微分方程式, plot(t,y) //绘出y对t的图形
35. function y=f(x),y=x*sin(30*x)/sqrt(1-((x/(2*%pi))^2)),endfunction //定义函数f -->intg(0,2*%pi,f) //对函数f积分,由0积分到2pi integrate('1/(1+(230*x-30)^2)','x',0,1) //對'1/(1+(230*x-30)^2)'作0到1之積分
36. apropos 文档中关键词搜寻;help 在线帮助;clear 从内存中清除变量和函数;exit 关闭SCILAB;quit 退出SCILAB; save 把内存变量存入磁盘, exec 运行脚本文件;mode 文件运行中的显示格式;getversion 显示SCILAB版本;ieee 浮点运算溢出显示模式选择;edit 文件编辑器;type 变量类型;what 列出SCILAB基本命令;format 设置数据输出格式;chdir 改变当前工作目录;getenv 给出环境值;mkdir 创建目录;pwd 显示当前工作目录;evstr 执行表达式;
37. 运算符,和特殊算符:*矩阵乘,.*数组乘,^ 矩阵乘方,.^ 数组乘方,/ 斜杠或右除,./或. 数组除,~= 不等号,~,not 逻辑非,' 复数转置号,.' 转置号,%eps 浮点误差容限, %i 虚数单位 = √(-1),%inf 正无穷大;%pi 圆周率, π=3.1415926535897....;
38. 基本数学函数: acos 反余弦;acosh 反双曲余弦;acot 反余切;exp 指数;log 自然对数;log10 常用对数;log2 以2为底的对数;abs 绝对值;conj 复数共轭;imag 复数虚部;real 复数实部;fix 向零方向取整;ceil 向上(正无穷大方向)取整;floor 向下(负无穷大方向)取整;round 四舍五入取整;gsort 降次排序;gamma gamma函数;interp 插值函数;interpln 线性插值函数;intsplin 样条插值函数;smooth 样条平滑函数;spline 样条函数;quarewave 方波函数;
39. 基本矩阵函数和操作: eye 单位阵;zeros 全零矩阵;ones 全1 矩阵;rand 均匀分布随机阵;genmarkov 生成随机Markov矩阵;linspace 线性等分向量;logspace 对数等分向量;logm 矩阵对数运算;cumprod 矩阵元素累计乘;cumsum 矩阵元素累计和;toeplitz Toeplitz 矩阵;disp 显示矩阵和文字内容;length 确定向量的长度;size 确定矩阵的维数;diag 创建对角阵或抽取对角向量;find 找出非零元素1的下标;matrix 矩阵变维;rot90 矩阵逆时针旋转90度;sub2ind 据全下标换算出单下标;tril 抽取下三角阵;triu 抽取上三角阵;conj 共轭矩阵;companion 伴随矩阵;det 行列式的值;norm 矩阵或向量范数;nnz 矩阵中非零元素个数;null 清空向量或矩阵中的某个元素;orth 正交基;rank 矩阵秩;trace 矩阵迹;cond 矩阵条件数;rcond 逆矩阵条件数;inv 矩阵的逆;lu LU分解或高斯消元法;pinv 伪逆;qr QR分解;givens Givens变换;linsolve 求解线性方程;lyap Lyapunov方程;hess Hessenberg 矩阵;poly 特征多项式;schur Schur 分解;expm 矩阵指数;expm1 矩阵指数的Pade逼近;expm2 用泰勒级数求矩阵指数;expm3 通过特征值和特征向量求矩阵指数;funm 计算一般矩阵函数;logm 矩阵对数 ;sqrtm 矩阵平方根 ;
40. 利用字符串实现公式的赋值计算:
调用格式:A=evstr(str)
描述: 该函数实现将原字符串str中的各个已知变量的数值带入公式进行计算。
例子:-->A=[‘x+y’ ‘y^w’;’z^2+2*z-3’ ‘w*v’]
A =
!x+y y^w !
! !
!z^2+2*z-3 w*v !
-->x=1;y=2;z=3;w=4;v=5;
-->evstr(A)
ans =
3. 16.
12. 20.