please can someone look at this code an find why my runge kutta code don't work at i expect(the coefficient are constant althought they schouldn't). )

if true
  % code

function yprime = gsm2(t,y) drehz = y(1); strom = y(2); global k1 n1 nn Mn Jm Ra Re Oh p Jg D La

 k1= 720;
 n1=50;
  Ui=[4.8 105 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 420 440]; %Induzierte Ankerspannung
  Ie= [0 0.2 0.245 0.29 0.34 0.375 0.415 0.46 0.505 0.555 0.6 0.65 0.705 0.76 0.84 0.905 1 1.115 1.27]; %Erregerstrom

%Hauptpolwindungszahl
We= 952*2*p;
%--------------------------------------------------
%Fluss
%--------------------------------------------------
c= 1/(k1*n1);
PHI= c.*Ui;  %Berechnung des Flusses %   Phi= [0.0001 0.0029 0.0033 0.0039 0.0044 0.005 0.0056 0.0061 0.0067 0.0072 0.0078 0.0083 0.0089 0.0094 0.01 0.0106 0.0111 0.0117 0.0122];

ie=0.36;
phi=spline(Ie,PHI,ie);  % PHI(Ie)-Kennlinie vom Fluss

%--------------------------------------------------
%Ableitung des Flusses
%--------------------------------------------------
i=1:18;
   Phi= PHI(i+1)-PHI(i);
   dIe= Ie(i+1)-Ie(i);
   k= Phi./dIe;
   PHi= [k 0.0035]; %Zuweisen einer zusätlichen Variable, sodass es kein Dimensionproblem bei der Benutzung von der Spline-Fkt gibt.

ie=0.36;
phI=spline(Ie,PHi,ie); % PHi(Ie)-Kennlinie der Ableitung des Flusses

%--------------------------------------------------
%DGL.-System 1. Ordnung des DC-Generators
%--------------------------------------------------

   a = n1-y(1);
   a1 = n1-nn;
   a2 = k1*phi*y(2);
   a3 = 2*pi*(Jg+Jm);
   a4= (a./a1)*Mn;
   a5= a2/(2*pi);
   dndt = (a4-a5)/a3;     %   drehzprime = ((((n1-y(1,:))/(n1-nn))*Mn)-((k1*phi*y(2,:))/(2*pi)))/(2*pi*(Jg+Jm));

   b= k1*phi*y(1);
   b1= (Ra+Re)*y(2);
   b2=(1+Oh)*We*phI;
   b3= La+b2;
   didt= (b-b1)/b3;     %   stromprime = ((k1*phi*y(1,:))-((Ra+Re)*y(2,:)))/(La+((1+Oh)*We*PHi));

yprime = [dndt; didt]; end

and the main to solve it is in the file. thx a lot

Answer 1

[kamdam Sophia]

sorry i have forgoteen to put the main clc global k1 n1 nn Mn Jm Ra Re Oh p Jg D La %-------------------------------------------------- % DC Motor Daten %-------------------------------------------------- D= 0.158; %Ankerdurchmesser p=2; %Polpaarzahl der Gsm l=0.25; %Blechpaketlänge rho=7874; %Dichte des Eisens k1= 720; %Maschinenkonstante Ra= 0.198; %Ankerkreiswiderstand La= 0.0042; %Ankerkreisinduktivität Re= 115; %Erregerwicklungswiderstand Oh=0.1; %Streuziffer

%-------------------------------------------------- % ASM Daten %-------------------------------------------------- n1=50; %Synchrone Drehzal nn= 48.8167; %Nenndrehzahl Mn=24; %Nenndrehmoment Jm= 0.024; %Massenträgheitsmoment der ASM

%-------------------------------------------------- %Berechnung der Massenträgheitsmoment(Jg) der GSM %-------------------------------------------------- m= rho*pi*((D*D)/4)*l; %Berechnung der Masse Jg= 0.5*m*((D*D)/4); %Berechnung der Massenträgheitsmoment der GSM

 %--------------------------------------------------
   %Runge-Kutta Lösungsalgorithmus
 %--------------------------------------------------
% Initial set up
h=2e-3;       % Step size of solution
t= 0:h:3*h;
y01=50;
y02=0;
y=[y01;y02];
t0= 0;
drehz0=50;        % Initial Drehz at t0
strom0=0;        % Initial strom at t0
% length(t)
% Create variables to store successive approximations
time  = t0;                 % stores all Time values
drehz = drehz0;             % stores all Drehz values
strom = strom0;             % stores all strom values

for n = 1:length(t)

   % Wrap Drehz and strom in y for convenience
   Y = @(t,x)gsm1(t,y(:,n));

   % Rechnung des k_Wertes für die nächste R-K Annäherung
    k_1n = Y(t(n),y(:,n));

  k_2n = Y(t(n)+0.5*h, y(:,n)+0.5*h*k_1n);

  k_3n = Y(t(n)+0.5*h, y(:,n)+0.5*h*k_2n);

  k_4n = Y(t(n)+h, y(:,n)+k_3n*h);

   % Rechnung der Annäherung für jede nächste ODE
    y_next = y(:,n)+(h/6)*(k_1n+2*k_2n+2*k_3n+k_4n);

   % Annäherung von der drehz und dem strom
       drehz = [drehz y_next(1)];
       strom = [strom y_next(2)];
       time  = [time time(n)+h];
       y=[y y_next];

end % Plot drehz vs time and strom versus time figure; plot(time,strom,'g');