run time for runge kutta with adaptive step size for System equations

Question

uzi le 22 Oct 2014

0
Lien

Utiliser le lien direct vers cette question

https://fr.mathworks.com/matlabcentral/answers/159695-run-time-for-runge-kutta-with-adaptive-step-size-for-system-equations

Réponse apportée : uzi le 24 Oct 2014

calc1.m

hi :-) i need to solve the system dy1/dx=f1(x,y1,...yn) dy2/dx=f2(x,y1,...yn) dy3/dx=f3(x,y1,...yn) ... ... dyn/dx=fn(x,y1,...yn) when n~~5000 so i wrote a runge kutta with adaptive step size

the equations come from "temp4000.txt" file and in "read_input" i convert the txt file to cell array. In "runge_kutta" i did the main calc.

my Problem is in the loops in function "calc_step", every loop Continue 4s!!

note: f is cell array for function f1,f2,f3,f4,f5,....,fn x is a number, y is a Vector y=[y1 y2 y3 y4....yn] when n~~5000 Thanks!!!

%%%%%%%%%%%%% %%%main function%% %%%%%%%%%%%%% function u=calc1(n,m,w0); x=read_input(); u=runge_kutta(x,n,m,w0); end

%%%%%%%%%%%%% %calc runge kutta%% %%%%%%%%%%%%% function Vector=runge_kutta(f,n1,n2,w0); % error for one step relative_error=0.01;

    n=size(f,2);
    %  if the num of Equations not EqualInitial conditions
    if n~=size(w0,2)        
        n
        size(w0,2)
        return
    end
    % set Initial data
    h=0.05;
    x=n1;
    y=w0;
    Vector=[x y];
    while x<n2
        x=x+h;
        y0=calc_step(f,x,y,h);                             % clac step for h
        y1=calc_step(f,x,y,h/2);                        % clac step for h/2
       error_step=calc_error_step(y0,y1);    % clac the error for the step
       OK=1;
       if error_step>relative_error
            i=n;
            x=x-h;
            h=0.9*h*(relative_error/error_step)^0.25;
            OK=0;
            x;
       end
       % Insert the  Solution to matrix
        if OK==1
            OK=1;
            y=y0;
            Vector=[Vector ; x y];
            h=0.9*h*(relative_error/error_step)^0.2;
            x
        end
    end
    Vector;
end

%%%%%%%%%%%%%%%%% % calc_step by runge kutta% %%%%%%%%%%%%%%%%% function y=calc_step(f,x,y,h); n=size(f,2); tic f1=zeros(1,n); g= @(l) l(x,y); f1=cellfun(g,f); toc % calc k2 tic x2=x+h/2; y2=y+h*f1/2; f2=zeros(1,n); parfor i=1:n f2(i)=f{i}(x2,y2); end toc % calc k3 tic x3=x+h/2; y3=y+h*f2/2; f3=zeros(1,n); parfor i=1:n f3(i)=f{i}(x3,y3); end toc % calc k4 tic x4=x+h/2; y4=y+h*f3; f4=zeros(1,n); parfor i=1:n f4(i)=f{i}(x4,y4); end toc y=y+h*(f1+2*f2+2*f3+f4)/6; end

%%%%%%%%%%%%%%% %calc relative error step% %%%%%%%%%%%%%%% function error=calc_error_step(y0,y1); error=abs((y0-y1)/y0); error=max(error); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %read the text file and return array (f_1 f_2 f_3) %the command to activate the function is array{i}(x,y) when y=[y_1 y_2.... y_n] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function input=read_input();

    %read the file
    input1=importdata('temp4000.txt');
    n=size(input1,1);
    %Conversion to string
    input2=cell2mat(input1(1));
    input3{1}=input2;
    for i=2:n
        input2=cell2mat(input1(i));
        input3{i}=input2;
    end
    %Creation a string '=@(x,y_1,...,y_n)'
    per=[];
    per=['=@(x,y)'];
    %add the string '=@(x,y_1,...,y_n)' to input3
    for i=1:n
        input3{i}=['f_' num2str(i) per input3{i}];
    end
    %activate the function
    for i=1:n
        eval(input3{i});
    end
    temp='input={';
    for i=1:n
        temp=[temp ' f_' num2str(i)];
    end
    temp=[temp '};'];
    eval(temp);
end

2 commentaires
Afficher AucuneMasquer Aucune

uzi le 22 Oct 2014

the "temp4000.txt"

Geoff Hayes le 22 Oct 2014

Since you have attached the calc1.m file, then you can remove the code from your question as it just clutters it up.

You may want to include the text file, temp4000.txt, so that someone else can try to run your code. As well, include the inputs to the calc1 function. What is n, m, and w0?

Have you tried to profile your code? See profile for details.

Connectez-vous pour commenter.

Connectez-vous pour répondre à cette question.