为了模拟，自学的matlab，心力憔悴，救救孩子。
原方程、文献图如图所示。其中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]); %求数值解
plot(xi,u(:,1),'-r',xi,u(:,2),'.k',xi,u(:,3),'*g',xi,u(:,4),':b')
积分方程组是基于上述结果的，也不知怎么搞：
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')
不知怎么发不了图片，图片：https://zhidao.baidu.com/questio ... ?entry=qb_uhome_tag
返回小木虫查看更多
京公网安备 11010802022153号
我计算了一下那几个N，算出是固定的一个数
你把代码发出来，我学习一下。我的还是运行不出来：
函数文件：
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)]; %等号左边的；等同等号右边的
end
end
---------
N = quad(@polymer,0,10);
R=1;
V=1;
a=0.1;
[xi,u]=ode45(@monomer,[0 10],[100;0;0;0]); %求数值解[100;0;0;0]为四条曲线的纵坐标初始值，等同于[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'); %通过图例'N_1','N_2','N_3','N_4'修改
xlabel('\xi');
ylabel('Amplitude')
-------------
运行结果：
>> polymer_s
索引超出数组元素的数目(0)。
出错 polymer (line 5)
N(1) = u(1)+1/2.*V.*u(2);
出错 quad (line 67)
y = f(x, varargin{:});
出错 polymer_s (line 1)
N = quad(@polymer,0,10);
>>
，
你私发整篇文献过来学习一下吧
我用cumtrapz解出来了
可以发过来学习一下吗？
发过去了