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

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

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

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

账号		自动登录	找回密码
密码			我要加入

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

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。