声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3065|回复: 0

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

[复制链接]
发表于 2006-5-9 15:57 | 显示全部楼层 |阅读模式

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

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

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

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-14 13:40 , Processed in 0.061819 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表