matlab调用龙格库塔法出错In an assignment A(I) = B,the number of elements in B and I must be the使用RK4求解常微分出错,常微方程如下,h+d=23为常量,欧米茄是已知的数组,共27个元素,程序中用omg定义.θ是定义的

matlab调用龙格库塔法出错In an assignment A(I) = B,the number of elements in B and I must be the使用RK4求解常微分出错,常微方程如下,h+d=23为常量,欧米茄是已知的数组,共27个元素,程序中用omg定义.θ是定义的
matlab调用龙格库塔法出错In an assignment A(I) = B,the number of elements in B and I must be the
使用RK4求解常微分出错,常微方程如下,h+d=23为常量,欧米茄是已知的数组,共27个元素,程序中用omg定义.θ是定义的等差数组,用seita定义,是自变量,应该也是27个元素,x是应变量.θ初始为0,x为-15 .


function inner = sax(varargin)
clc,clear
h1=8;d=15;l=2000;n=1.4;
h=pi/160;max_omg=30*pi/180;xn=pi/6;y0=0;
seita=[0:pi/160:pi/6];
lk=length(seita);
yt=[0,45.3421,90.6667,135.9563,181.1935,226.3609,271.441,316.4164,361.2699,405.9841,450.5418,494.9258,539.119,583.1043,626.8649,670.3837,713.6442,756.6295,799.3231,841.7086,883.7695,925.4898,966.8532,1007.844,1048.4461,1088.6441,1128.4224]
omg=[0,0.0152,0.0303,0.0454,0.0605,0.0754,0.0902,0.1048,0.1193,0.1336,0.1477,0.1615,0.1752,0.1885,0.2016,0.2144,0.227,0.2392,0.2511,0.262,0.274,0.285,0.2957,0.3061,0.3161,0.3258,0.3353]
global omg;global seita;


[y,x]=lgkt(d,lk,h) ;%调用rk4


function z=f(seita,x) %rk4子函数
global omg;global seita;
h1=8;d=15;n=1.4;
fenzi=h1+d+x;
fenmu=(-(n*cos(omg)-cos(seita))/(n*sin(omg)-sin(seita))-tan(seita)).*((cos(seita)).^2);
z=fenzi./fenmu;%常微分方程
function [seita,x]=lgkt(d,lk,h)
global seita;
x(1)=-d;
y1=seita;
y1(1)=0;
for i=1:(lk-1)
    K1=f(y1(i),x(i));     %x为x左边 y为出射角度θ
    K2=f(y1(i)+h/2,x(i)+h/2*K1);
    K3=f(y1(i)+h/2,x(i)+h/2*K2);
    K4=f(y1(i)+h,x(i)+h*K3);
    x(i+1)=x(i)+h/6*(K1+2*K2+2*K3+K4);
end
y=y1;
plot(x,y)


matlab提示
  In an assignment  A(I) = B, the number of elements in B and
 I must be the same.
Error in ==> axs>lgkt at 31
    x(i+1)=x(i)+h/6*(K1+2*K2+2*K3+K4);
Error in ==> axs at 13
[y,x]=lgkt(d,lk,h) ;

matlab调用龙格库塔法出错In an assignment A(I) = B,the number of elements in B and I must be the使用RK4求解常微分出错,常微方程如下,h+d=23为常量,欧米茄是已知的数组,共27个元素,程序中用omg定义.θ是定义的
我只能说你的程序很乱,一堆没有用到的变量,不知道你在想什么.最关键的是常微分方程的omg不应该是定值吗,也即是说对每一个给定的omg值都有一个常微分方程的解与之对应.我帮你改了一下.仅供参考.
function sax
clc;clear;
h=pi/160;
seita=0:h:pi/6;
x0=-15;
x=lgkt(seita,x0,h);%调用rk4
plot(seita,x)
end
function x=lgkt(seita,x0,h)
y1=seita;
x=zeros(size(y1));
x(1)=x0;
for i=1:length(seita)-1
K1=f(y1(i),x(i));%x为x左边y为出射角度θ
K2=f(y1(i)+h/2,x(i)+h/2*K1);
K3=f(y1(i)+h/2,x(i)+h/2*K2);
K4=f(y1(i)+h,x(i)+h*K3);
x(i+1)=x(i)+h/6*(K1+2*K2+2*K3+K4);
end
end
function z=f(seita,x)%rk4子函数
h1=8;%这里修改
d=15;
n=1.4;
omg=0.2;
fenzi=h1+d+x;
fenmu=(-(n*cos(omg)-cos(seita))/(n*sin(omg)-sin(seita))-tan(seita)).*((cos(seita)).^2);
z=fenzi./fenmu;%常微分方程
end