科学分析软件 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

grass22 发表于 2006-5-9 15:57

MATHEMATICA讲座(2) ----徐安农教授

MATHEMATICA讲座第六讲 线性方程组的表达方式和解法 一.向量和矩阵的输入 1.Range Range Range 执行算例 Range Range 2.Table[(1/2)^n,{n,0,10}] (*由通项构造表*) Table[{f1(n),F2(n)},{n,n1,n2,h}] 执行算例 Table[(1/2)^n,{n,0,10}] A=Table,{x,0,n}]],{n,1,11,2}] (*Sin的六个正规化后的幂级数展开式的表*) 执行算例 Table[{x,x^2},{x,-1.0,1.0,0.2}] Table,{i,2,5}] Table,{n,0,17}] Table,{n,1,5}] Table,{10}] Table,{10}] 3.Array[函数,n] (*由函数表达式构造表*) Array[函数，{n1,n2,n3,...}] 执行算例 Array（*等价于Table,{x,5}]*） Array （*等价于Table{n,0,10},{m,0,10}]*） 4.NestList,n0,k] 用递推公式建立表元素 Clear t1=(Sqrt+1)/2 t2=(1-Sqrt)/2 fnt=Table[(t1^(n+1)-t2^(n+1))/Sqrt,{n,0,40}]//N rnt=Table]/fnt[],{n,2,12}] 5.向量与矩阵的标准输入法 A={x1,x2,x3,...} A={{a11,a12,a13},{a21,a22,a23},{a31,a32,a33}} ColumnForm[{a1,a2,...,an}]把向量用列方式输出 MatrixForm 用矩阵方式显示 IdentityMatrix 生成n阶单位阵 DiagonalMatrix[{a11,a22,...,ann}]生成对角阵 执行算例 A1={1,2,3} ColumnForm A2=DiagonalMatrix[{1,2,3,4,5}] MatrixForm IdentityMatrix DiagonalMatrix[{1,2,3,4,5}] MatrixForm[%] 二.行列式 Det 执行算例 A={{1,2,3},{4,5,6},{7,8,9}} MatrixForm Det 三.矩阵求逆,求特征值，特征向量 Inverse Eigenvalue Eigenvector 执行算例 Clear A={{4,6,0},{-3,-5,0},{-3,-6,1}} Eigenvalues Eigenvectors 四.恰定方程求解 问题1 x1+6x2+36x3=104 x1+10x2+100x3=160 x1+20x2+400x3=370 程序 Clear A1={{1,6,36},{1,10,100},{1,20,400}}; b={104,160,370}; LinearSolve(*求方程组的解*) X=Inverse.b (*用求逆矩阵方法求解*) 五.欠定方程求解 问题2 2x1+x2-x3+x4=1 x1+2x2+x3-x4=2 x1+x2+2x3+x4=3 程序 A2={{2,1,-1,1},{1,2,1,-1},{1,1,2,1}} b2={1,2,3}; X={x1,x2,x3,x4}; Solve MATHEMATICA讲座第七讲 函数的插值 一.拉格朗日插值 L={List} InterpolatingPolynomial 执行算例1 两点线性插值 L={{0,0.3},{0.2,0.45}} I=InterpolatingPolynomial 执行算例2 三点抛物插值 L1={{0,0.3},{0.2,0.45},{0.4,0.15}} I1=InterpolatingPolynomial 执行算例3 多点拉格朗日插值 L2={{0,0.3},{0.2,0.45},{0.3,0.47}, {0.52,0.50},{0.64,0.38},{0.7,0.33},{1.0,0.24}} I2=InterpolatingPolynomial Plot[%,{x,-0.25,1.05}] 执行算例4 作正弦在0，P上五点插值函数图形 g0=Plot,{x,0,Pi}] L=Line},{x,0,Pi,Pi/4}]] g=Graphics Show sinAp:=Graphics[{Line}, {x,0,Pi,Pi/(n+1)}]]}] sinAp Show 二.龙格现像演示 L=Table[{x,1/(1+25*x^2)},{x,-1,1,0.2}] a=InterpolatingPolynomial b=Plot}] c=Plot Show 三. 两点三次Hermite插值 执行算例5 Clear x1={0,1};y1={1,2};m={1/2,1/2}; h0:=(1+2*x)*(x-1)^2; h1=(1-2(x-1))*(x/(x-1))^2 H0=x*(x-1)^2 H1=(x-1)*x^2 H=y1[]*h0+y1[]*h1 +m[]*H0+m[]*H1 %/.{x->0.55} 四. N+1个节点的2N+1次Hermite插值 执行算例6 Clear x0={0.4,0.5,0.6,0.7,0.8} y=Table,{x,0.40,0.80,0.10}] m=Table bb=InterpolatingPolynomial Simplify] bb w=(x-x0[])*(x-x0[])*(x-x0[])*(x-x0[])*(x-x0[]) w1=D,x]; Simplify] w2=D,{x,2}]; Simplify] For[ i=1,i<=5,i++, L:=w/((x-x0[])*w1]])]; h:=(1-w2]]*(x-x0[])/w1]])*L^2; H:=L^2*(x-x0[]); ] Hm=Sum]*h+m[]*H,{i,1,5,1}]; Simplify] Hm 拟合 一.一元线性拟合 执行算例 b2={{100,45},{110,51},{120,54},{130,61},{140,66}, {150,70},{160,74},{170,78},{180,85},{190,89}} fp=ListPlot, RGBColor}] ft1=Fit gp=Plot}]; Show 二.抛物线拟合 执行算例 B=Table,{n,20}] t1=ListPlot, PointSize}] f=Fit t2=Plot Show 三.多项式拟合 执行算例 data={{0,1.2},{1,1.4},{2,1.3},{3,1.5},{4,1.3},{5,1.3},{6,1.1}}; t2=ListPlot, RGBColor}] fx=Fit t1=Plot] Show 执行算例 b3={{1,4},{2,6.4},{3,8.0},{4,8.4},{5,9.28},{6,9.5}, {7,9.7},{8,9.86},{9,10.0},{10,10.2},{11,10.32}, {12,10.42},{13,10.5},{14,10.55},{15,10.58}, {16,10.6}} gp=ListPlot, PointSize}] ft2=Fit,x] fp=Plot}] Show 四.非线性拟合---指数拟合 执行算例求一个经验函数,型如 x 1 2 3 4 5 6 7 8 y 15.3 20.5 27.4 36.5 49.1 65.5 87.8 117.6 程序 b4={{1,15.3},{2,20.5},{3,27.4},{4,36.6},{5,49.1}, {6,65.5},{7,87.8},{8,117.6}}; gb4=ListPlot, PointSize}] y4=Table],{i,1,8}]; ly4=Log; fy4=Fit s:=Exp ty=Plot,{x,1,8},Axes->{2,60}, AspectRatio->1, PlotStyle->{RGBColor}, PlotRange->{10,120}] Show 执行算例求一个经验函数,型如y=a*exp(-bx)与所给数据拟合 x 0.4 0.5 0.6 0.7 y 1.75 1.34 1.00 0.74 程序 Clear fx:=x fy:=Log biao={{0.4,1.75},{0.5,1.34},{0.6,1.00},{0.7,0.74}}; nb=Table[{fx]],fy]]},{i,1,4}]; (*拟合方程*) ft=Fit ft1=Exp (*拟合曲线*) t1=Plot}] t2=ListPlot,PointSize}] Show 五.用正交多项式作拟合 区间上的勒让得多项式 (*定义勒让得函数(n=10)*) Clear n=10 P0:=1 P1:=1-2x/n P2:=1-6 x/n+6 x*(x-1)/(n*(n-1)) P3:=1-12 x/n+30x*(x-1)/(n*(n-1))-20 x*(x-1)*(x-2)/(n*(n-1)*(n-2)) P4:=1-2 x+x*(x-1)-140x*(x-1)*(x-2)/720+70x*(x-1)*(x-2)*(x-3)/(10*9*8*7) P5:=1-30x/n+210x*(x-1)/(n*(n-1))-560x*(x-1)*(x-2)/(n*(n-1)*(n-2)) +630x*(x-1)*(x-2)*(x-3)/(n*(n-1)*(n-2)*(n-3)) -252 x*(x-1)*(x-2)*(x-3)*(x-4)/(n*(n-1)*(n-2)*(n-3)*(n-4)) (*输入初始数据*) t={0,5,10,15,20,25,30,35,40,45,50}; y={0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.60,4.66}; (*做变量替换*) x=t/5; (*计算各多项式在节点处的值*) A={P0,P1,P2,P3,P4,P5} (*计算每一行元素平方的和*) s=Table; For]=0; For ]=s[]+A[]^2] ] N (*计算Pk(xi)*yi*) r=Table For]=b0[]*y[] ] For]=b0[]*y[] ] For]=b0[]*y[] ] For]=b0[]*y[] ] For]=b1[]*y[] ] For]=b2[]*y[] ] For]=b3[]*y[] ] For]=b4[]*y[] ] For]=b5[]*y[] ] r MATHEMATICA讲座第八讲 线性规划与非线性规划 模型 min(or max) f=c'X=c1x1+c2x2+...+cnxn s.t. MX>=b,X>=0 命令格式 ConstrainedMax ConstrainedMin LinearProgramming 执行算例1 求解线性规划问题 max f=2x+3y s.t. x+2y<=8 0<=x<=4,0<=y<=3 ConstrainedMax 第一个参数是目标函数的最优值，第二个参数是决策变量的取值。 执行算例2 min f=x+2y+3z s.t. 2x-y =1 x +z=1 x,y,z>=0 (*解法一*) ConstrainedMin (*解法二*) c={2,-3} M={{-1,-1},{1,-1},{1,0}} b={-10,2,1} q=LinearProgramming 执行算例3 求解线性规划问题 min -f =- 0.40x1-0.28x2-0.32x3-0.72x4-0.64x5-0.61x6 s.t. -0.01x1-0.01x2-0.01x3-0.03x4-0.03x5-0.03x6>=-850 -0.02x1 -0.05x4 >=-700 -0.02x2 -0.05x5 >=-100 -0.03x3 -0.08x6>=-900 x1,x2,x3,x4,x5,x6>=0 解： c={-0.4,-0.28,-0.32,-0.72,-0.64,-0.6}; A={{-0.01,-0.01,-0.01,-0.03,-0.03,-0.03}, {-0.02,0,0,-0.05,0,0},{0,-0.02,0,0,-0.05,0}, {0,0,-0.03,0,0,-0.08}}; b={-850,-700,-100,-900}; LinearProgramming (* 此命令只给出决策变量。*) 线性约束条件下的非线性规划问题 线性逼近法(FW法) 模型: (NLP)minf(x) S.t.: E={x|AX>=b,X>=0} 执行算例4 求解非线性规划问题 用线性逼近法求解非线性规划问题 目标函数 f(x1,x2)=(x1-1)^2+(x2-2)^2 约束条件 0<=x1<=2,0<=x2<=3 下面是第一次迭代 Clear f:=(x1-1)^2+(x2-2)^2; gradf={D,x1],D,x2]}; c={0.7,1.25}; (*C即基可行解X0*) s=gradf/.{x1->Part,x2->Part} (*求X0点梯度*) {-0.6, -1.5} x1=.; (*清除X1,X2的值*) x2=.; p:=s.{u,v} a=ConstrainedMin,{x1<=2,x2<=3},{x1,x2}] (*求出最优解Y0*) b={x1,x2}/.Part; e=b-c; pf=s.e; If>0.01, g=c+d*e; t:=f,Part]; w=FindMinimum,{d,1,0,1}]; c=c+(d/.Part)e; ] Print["c1=",c] (* 得到初始点c1,将其替换c,运算后的结果继续代替C,直到w的绝对 值小于0.01为止。*) *)

页: [1]

声振论坛's Archiver

MATHEMATICA讲座(2) ----徐安农教授