為了模擬，自學(xué)的matlab，心力憔悴，救救孩子。
原方程、文獻(xiàn)圖如圖所示。其中R=1;V=1;a=0.1;我模擬的微分模擬圖不知道為啥是這樣的。
微分方程組：
function du=f(xi,u)
%a替換alpha
syms R V xi a u0 u1 u2 u3
du=zeros(4,1);
du0=-u0-1./2*V*u1; %22 初始值
du1=-V*u0-(1+R*exp(-a*xi)*cosh(a*V *xi))*u1-(1./2*V-R*exp(-a*xi)*sinh(a*V *xi))*u2; %23;a替換alpha
du2=-(1./2*V-R*exp(-a*xi)*sinh(a*V *xi))*u1-(1+4*R*exp(-a*xi)*cosh(a*V *xi))*u2-(1./2*V-3*R*exp(-a*xi)*sinh(a*V *xi))*u3; %24
du3=-(1./2*V-R*exp(-a*xi)*sinh(a*V *xi))*u2-(1+9*R*exp(-a*xi)*cosh(a*V *xi))*u3'; %25
end
R=1;
V=1;
a=0.1;
[xi,u]=ode45('f',[0 10],[100 0 0 0]); %求數(shù)值解
plot(xi,u(:,1),'-r',xi,u(:,2),'.k',xi,u(:,3),'*g',xi,u(:,4),':b')
積分方程組是基于上述結(jié)果的，也不知怎么搞：
function N=integral(xi)
syms xi V N u0 u1 u2 u3 N0 N1 N2 N3
N0=u0+1./2*V*u1;
N1=V*u0+u1+1./2*V*u2;
N2=1./2*V*u1+u2+1./2*V*u3;
N3=1./2*V*u2+u3
end
N0=quad(@integral,0,10);
N1=quad(@integral,0,10);
N2=quad(@integral,0,10);
N3=quad(@integral,0,10);
plot(xi,N(:,1),'-r',xi,N(:,2),'.k',xi,N(:,3),'*g',xi,N(:,4),':b')
不知怎么發(fā)不了圖片，圖片：https://zhidao.baidu.com/questio ... ?entry=qb_uhome_tag
返回小木蟲(chóng)查看更多
京公網(wǎng)安備 11010802022153號(hào)
我計(jì)算了一下那幾個(gè)N，算出是固定的一個(gè)數(shù)
你把代碼發(fā)出來(lái)，我學(xué)習(xí)一下。我的還是運(yùn)行不出來(lái)：
函數(shù)文件：
function N = polymer(xi)
global V u;
N(1) = u(1)+1/2.*V.*u(2);
N(2) = V.*u(1)+u(2)+1/2.*V.*u(3);
N(3) = 1/2.*V.*u(2)+u(3)+1/2.*V.*u(4);
N(4) = 1/2.*V.*u(3)+u(4);
N=@monomer;
function u=monomer(xi,u)
%a替換alpha
global R a;
du(1)=-u(1)-1./2*V*u(2); %22 初始值
du(2)=-V*u(1)-(1+R*exp(-a*xi)*cosh(a*V *xi))*u(2)-(1./2*V-R*exp(-a*xi)*sinh(a*V *xi))*u(3); %23;a替換alpha
du(3)=-(1./2*V-R*exp(-a*xi)*sinh(a*V *xi))*u(2)-(1+4*R*exp(-a*xi)*cosh(a*V *xi))*u(3)-(1./2*V-3*R*exp(-a*xi)*sinh(a*V *xi))*u(4); %24
du(4)=-(1./2*V-R*exp(-a*xi)*sinh(a*V *xi))*u(3)-(1+9*R*exp(-a*xi)*cosh(a*V *xi))*u(4); %25
u=[du(1);du(2);du(3);du(4)]; %等號(hào)左邊的；等同等號(hào)右邊的
end
end
---------
N = quad(@polymer,0,10);
R=1;
V=1;
a=0.1;
[xi,u]=ode45(@monomer,[0 10],[100;0;0;0]); %求數(shù)值解[100;0;0;0]為四條曲線的縱坐標(biāo)初始值，等同于[100 0 0 0]
u=polymer(xi,u);
plot(xi,N(:,1),'-k',xi,N(:,2),':c',xi,N(:,3),'--g',xi,N(:,4),'-.r','LineWidth',2);
legend('N_0','N_1','N_2','N_3'); %通過(guò)圖例'N_1','N_2','N_3','N_4'修改
xlabel('\xi');
ylabel('Amplitude')
-------------
運(yùn)行結(jié)果：
>> polymer_s
索引超出數(shù)組元素的數(shù)目(0)。
出錯(cuò) polymer (line 5)
N(1) = u(1)+1/2.*V.*u(2);
出錯(cuò) quad (line 67)
y = f(x, varargin{:});
出錯(cuò) polymer_s (line 1)
N = quad(@polymer,0,10);
>>
，
你私發(fā)整篇文獻(xiàn)過(guò)來(lái)學(xué)習(xí)一下吧
我用cumtrapz解出來(lái)了
可以發(fā)過(guò)來(lái)學(xué)習(xí)一下嗎？
發(fā)過(guò)去了