matlab软件编程物理模型

%理想d簧阵子简谐运动

%Clear

rectangle('position',[12,8.5,2,0.3],'FaceColor',[0.5,0.3,0.4])

axis([0,15,-1,10])

%画顶板

hold on

plot([13,13],[7,8.5],'r','linewidth',2)

%画直线

y=2:.2:7

M=length(y)

x=12+mod(1:M,2)*2

x(1)=13

x(end-3:end)=13

D=plot(x,y)

%d簧

C=0:.1:2*pir=0.35

t1=r*sin(C)

F1=fill(13+r*cos(C),2+t1,'r')

% 球

set(gca,'ytick',[0:2:9])

set(gca,'yticklabels',num2str([-1:3]'))

plot([0,15],[3.3,3.3],'black')

H1=plot([0,13],[3.3,3.3],'y')

% 句柄[黄线]

Q=plot(0,3.8,'color','r')

% 运动曲线

td=[]yd=[]

T=0

text(2,9,'理想中的d簧振子简谐振动','fontsize',16)

set(gcf,'doublebuffer','on')

while T<12

pause(0.2)

Dy=(3/2-1/2*sin(pi*T))*1/2

Y=-(y-2)*Dy+7

Yf=Y(end)+t1

td=[td,T]yd=[yd,Y(end)]

set(D,'ydata',Y)

set(F1,'ydata',Yf,'facecolor',rand(1,3))

set(H1,'xdata',[T,13],'ydata',[Y(end),Y(end)])

set(Q,'xdata',td,'ydata',yd)

T=T+0.1

end

具体仿真效果参见下面链接，这个是d簧阵子的动画效果，链接中还有其他集中模型的动画效果仿真

你好，首先我要说你选择matlab这一强大软件是明智的，它的功能十分全面，其优化工具箱解决你的问题十分方便线性规划的求解程序名为linprog，调用格式为[x,fval,exitflag,output,lambda] = linprog(c,A,b,Aeq,beq,lb,ub,x0,options)其中，x：最优解；val：最优解处的函数值；exitflag：程序结束时的状态指示（>0： 收敛，0： 函数调用次数或迭代次数达到最大值（该值在options中指定） <0： 不收敛）； Output：包含以下数据的结构变量（iterations 实际迭代次数， cgiterations 实际PCG迭代次数（大规模计算用），algorithm 实际使用的算法）；lambda：包含以下数据的结构变量（ineqlin 不等式约束的LAGRANGE乘子， eqlin 等式约束的LAGRANGE乘子， upper 上界约束的LAGRANGE乘子， lower 下界约束的LAGRANGE乘子）；c：目标函数矩阵；A/Aeq：不等式/等式限制条件系数矩阵；b/beq：不等式/等式限制条件常数项矩阵；lb：自变量定义域下限；ub：自变量定义域上限；x0：初始解（缺省时程序自动取x0=0）options：包含算法控制参数的结构 实际应用是并没有这么复杂，很多参数缺省就可以了，用你的例题演示就是：（在命令窗口输入）f=[-2-3]A=[0,14,21,1]b=[12206]lb=zeros(3,1)[x,fval]=linprog(f,A,b,[],[],lb)计算出x=[06]，即x1=0，x2=6fval=-18解释一下linprog函数只能求最小值，所以将目标函数系数全变为相反数，最终得到的结果应为fval的相反数希望以上内容对你学习matlab能有一定帮助，以后可以多看看帮助文件，里面讲解很详细。

1. 等量异号点电荷电场中电势分布

[x,y]=meshgrid(-2:0.1:2,-2:0.1:2)

%以0.1为步长建立平面数据网格

z=1./sqrt((x-1).^2+y.^2+0.01)... %写出电势表达式

-1./sqrt((x+1).^2+y.^2+0.01)

[px,py]=gradient(z)

%求电势在x,y方向的梯度即电场强度

contour(x,y,z,[-12,-8,-5,-3,-1,... %画出等势线

-0.5,-0.1,0.1,0.5,1,3,5,8,12])

hold on %作图控制

quiver(x,y,px,py,'k') %画出各点上电场的大小和方向

2. 等量同号点电荷的电场线的绘制

下面是写微分方程的函数文件：

function ydot=dcx1fun(t,y,flag,p1,p2)

%p1,p2是参量，表示电量

ydot=[p1*(y(1)+2)/(sqrt((y(1)+2).^2+y(2).^2).^3)+...

p2*(y(1)-2)/(sqrt((y(1)-2).^2+y(2).^2).^3)

%dx/dt=Ex

p1*y(2)/(sqrt((y(1)+2).^2+y(2).^2).^3)+...

p2*y(2)/(sqrt((y(1)-2).^2+y(2).^2).^3)]

%dy/dt=Ey

编写好函数文件后，命名为dcx1fun.m存在当前路径下，然后开始编写解微分方程的主程序dcx1.m：

p1=10p2=10%点电荷所带电量

axis([-5,5,-5,5])%设定坐标轴范围 -5≤x≤5,-5≤y≤5

hold on %图形控制,不可擦除模式

plot(2,0,'*r')plot(-2,0,'*r') %绘制两源电荷

a=(pi/24):pi/12:(2*pi-pi/24)

%圆周上电场线起点所对应的角度

b=0.1*cos(a)c=0.1*sin(a)

%电场线起点所对应的相对坐标

b1=-2+bb2=2+b%把起点圆周的圆心放置在源电荷处

b0=[b1 b2]c0=[c c]%初始条件，所有电场线的起点

%的横、纵坐标构成了矢量b0和c0

for i=1:48 %循环求解48次微分方程

[t,y]=ode45('dcx1fun',[0:0.05:40],[b0(i),c0(i)],[ ],p1,p2)

%调用ode45求解，对应一个初条件（起点），求解出一条电场线

plot(y(:,1),y(:,2),'b') %绘制出此条电场线

end %结束循环，共绘制出48条电场线

参考于《Matlab 在基础物理学中的应用》

http://bnucourse.bnu.edu.cn/course/physics/05/jcwlxshyjy.pdf

---