| 查看: 1099 | 回复: 6 | ||
[求助]
请大牛看看这个程序哪里出错了 很急 已有2人参与
|
| 这个用二分法解的方程 只能当f=50的时候能算出解 换成100 其他的数就是无解 有知道为什么的大神么 可以是有偿的 很急 很重要 谢谢 |
回复此楼» 猜你喜欢
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有23人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急
已经有15人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急
已经有6人回复
-
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急
已经有5人回复
-
售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急
已经有5人回复
-
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急
已经有6人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有6人回复
-
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急
已经有7人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有7人回复
-
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急
已经有10人回复
2楼2015-12-28 09:57:34
getengqing
木虫 (正式写手)
- 应助: 104 (高中生)
- 金币: 3759.8
- 散金: 30
- 红花: 5
- 帖子: 592
- 在线: 110.5小时
- 虫号: 2986702
- 注册: 2014-02-23
- 性别: GG
- 专业: 自然语言理解与机器翻译
3楼2015-12-28 10:22:49
4楼2016-01-02 19:08:13
5楼2016-01-02 19:24:35
|
%% clear; r0=5e-3; % r0=1e-3; % miu=8; % sigma_e=1/(60e-8); f=100;T=1/f; H0=40; w=2*pi*f; H=@(t)40*sin(2*pi*f*t); miu_0=4*pi*1e-7; d=10e-3; pho=60e-8; S=78.5e-6; beta=16; V_0=1e3; G=0.1356; k_e=3.283; alpha=-0.02; c=0.18; Ms=770; a=7.012; %% 欧拉法求数值解 I_max=200; dt=T/I_max;t0=0; %% 起点值; M=0; c1=miu_0*d.^2/(2*pho*beta); c2=sqrt(miu_0*G*S*V_0/(2*pho*beta)); t=t0; % H=@(x,y)H0*exp(j*w*t+sqrt(j*miu*sigma_e*w).*(sqrt(x.^2+y.^2)-r1)).*r1./sqrt(x.^2+y.^2); % ymax=@(x)sqrt(r1^2-x.^2); % ymin=@(x)sqrt(r0^2-x.^2); % HH=quad2d(H,0,r0,ymin,ymax)+quad2d(H,r0,r1,0,ymax); % HH=real(HH)*4/(pi*(r1^2-r0^2)); HH=0; for i=1:I_max; t=t0+i*dt; % H=@(x,y)H0*exp(j*w*t+sqrt(j*miu*sigma_e*w).*(sqrt(x.^2+y.^2)-r1)).*r1./sqrt(x.^2+y.^2); % Ht=quad2d(H,0,r0,ymin,ymax)+quad2d(H,r0,r1,0,ymax); % Ht=real(Ht)*4/(pi*(r1^2-r0^2)); Ht=H(t); dH=Ht-HH(i); %% mean_H=quad2d "积分" if Ht<0 delta=-1; else delta=1; end HH(i+1)=Ht; M(i+1)=erfen(c1,dH,dt,c2,k_e,delta,M(i),alpha,Ht); B(i+1)=miu_0*(M(i+1)+HH(i+1)); end plot(HH,B,'ro') xlabel('磁场强度') ylabel('磁感应强度') function M_1=erfen(c1,dH,dt,c2,k_e,delta,M,alpha,H) fun=@(x)odeff(x,c1,dH,dt,c2,k_e,delta,M,alpha,H); left=-1000; right=1000; n=1000;dh=(right-left)/n; for i=1:n; x_l=left+(i-1)*dh; x_r=left+i*dh; if fun(x_l)*fun(x_r)>0; continue;end if fun(x_l)==0;M_1=x_l;break;end; if fun(x_r)==0;M_1=x_r;break;end; d=1; while d>1e-5; x=(x_l+x_r)/2; if fun(x_l)*fun(x)==0; M_1=x;break; elseif fun(x_l)*fun(x)<0; x_r=x; else x_l=x; end d=abs(x_l-x_r); end M_1=(x_l+x_r)/2; end end function ff=odeff(M_1,c1,dH,dt,c2,k_e,delta,M,alpha,H) c=0.18; Ms=770; a=7.012; dM=M_1-M; He=H+alpha*(M+dM); He_0=H-dH+alpha*M; dHe=He-He_0; M_an=Ms*(coth(He/a)-a/He); dM_an=M_an-Ms*(coth(He_0/a)-a/He_0); ff=c1*dH/dt*(dM/dH).^2+c2*sqrt(dH/dt)*(dM/dH).^1.5+... (k_e*delta-alpha*(M_an-M_1+k_e*delta*c*dM_an/dHe))*dM/dH-... (M_an-M_1+k_e*delta*c*dM_an/dHe); end |
6楼2016-01-02 19:33:07
|
%% clear; r0=5e-3; % r0=1e-3; % miu=8; % sigma_e=1/(60e-8); f=100;T=1/f; H0=40; w=2*pi*f; H=@(t)40*sin(2*pi*f*t); miu_0=4*pi*1e-7; d=10e-3; pho=60e-8; S=78.5e-6; beta=16; V_0=1e3; G=0.1356; k_e=3.283; alpha=-0.02; c=0.18; Ms=770; a=7.012; %% 欧拉法求数值解 I_max=200; dt=T/I_max;t0=0; %% 起点值; M=0; c1=miu_0*d.^2/(2*pho*beta); c2=sqrt(miu_0*G*S*V_0/(2*pho*beta)); t=t0; % H=@(x,y)H0*exp(j*w*t+sqrt(j*miu*sigma_e*w).*(sqrt(x.^2+y.^2)-r1)).*r1./sqrt(x.^2+y.^2); % ymax=@(x)sqrt(r1^2-x.^2); % ymin=@(x)sqrt(r0^2-x.^2); % HH=quad2d(H,0,r0,ymin,ymax)+quad2d(H,r0,r1,0,ymax); % HH=real(HH)*4/(pi*(r1^2-r0^2)); HH=0; for i=1:I_max; t=t0+i*dt; % H=@(x,y)H0*exp(j*w*t+sqrt(j*miu*sigma_e*w).*(sqrt(x.^2+y.^2)-r1)).*r1./sqrt(x.^2+y.^2); % Ht=quad2d(H,0,r0,ymin,ymax)+quad2d(H,r0,r1,0,ymax); % Ht=real(Ht)*4/(pi*(r1^2-r0^2)); Ht=H(t); dH=Ht-HH(i); %% mean_H=quad2d "积分" if Ht<0 delta=-1; else delta=1; end HH(i+1)=Ht; M(i+1)=erfen(c1,dH,dt,c2,k_e,delta,M(i),alpha,Ht); B(i+1)=miu_0*(M(i+1)+HH(i+1)); end plot(HH,B,'ro') xlabel('磁场强度') ylabel('磁感应强度') function M_1=erfen(c1,dH,dt,c2,k_e,delta,M,alpha,H) fun=@(x)odeff(x,c1,dH,dt,c2,k_e,delta,M,alpha,H); left=-1000; right=1000; n=1000;dh=(right-left)/n; for i=1:n; x_l=left+(i-1)*dh; x_r=left+i*dh; if fun(x_l)*fun(x_r)>0; continue;end if fun(x_l)==0;M_1=x_l;break;end; if fun(x_r)==0;M_1=x_r;break;end; d=1; while d>1e-5; x=(x_l+x_r)/2; if fun(x_l)*fun(x)==0; M_1=x;break; elseif fun(x_l)*fun(x)<0; x_r=x; else x_l=x; end d=abs(x_l-x_r); end M_1=(x_l+x_r)/2; end end function ff=odeff(M_1,c1,dH,dt,c2,k_e,delta,M,alpha,H) c=0.18; Ms=770; a=7.012; dM=M_1-M; He=H+alpha*(M+dM); He_0=H-dH+alpha*M; dHe=He-He_0; M_an=Ms*(coth(He/a)-a/He); dM_an=M_an-Ms*(coth(He_0/a)-a/He_0); ff=c1*dH/dt*(dM/dH).^2+c2*sqrt(dH/dt)*(dM/dH).^1.5+... (k_e*delta-alpha*(M_an-M_1+k_e*delta*c*dM_an/dHe))*dM/dH-... (M_an-M_1+k_e*delta*c*dM_an/dHe); end |
7楼2016-01-02 19:33:52
