Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Problem with matlab or maple

  1. Jul 1, 2005 #1
    hello
    i want to solve the folowing problem with matlab or maple but the answer is wrong:
    we have:

    rho0:=0.57;
    eta:=rho0*(evalf(Pi))/6;
    lambda1:=((1+2*eta)^2)/(1-eta)^4;
    lambda2:=-((1+eta/2)^2)/(1-eta)^4;
    c:=x->-(evalf(Pi))*((lambda1*(1-x^2))+4*eta*lambda2*(1-x^3)+(lambda1*eta/5)*(1-x^5));
    rho:=x->rho0*int(c(x-x1)*rho(x1),x1=0..20)+rho0+((rho0)^2)*c(0);
    rho(x1):=rho0;

    I want to solve it using iteration method. I mean I must do this until
    (abs(sum('t','i'=1..i)))^2 < 10^(-4):

    rho(x1):=rho0;
    if ((abs(sum('rho(x)-rho(x1)))^2<10^(-4)) then rho(x) is true else rho(x1)=rho(x) and repeat again.

    maple code that i was wrote:

    Code (Text):
    rho0:=0.57;
    eta:=rho0*(evalf(Pi))/6;
    lambda1:=((1+2*eta)^2)/(1-eta)^4;
    lambda2:=-((1+eta/2)^2)/(1-eta)^4;
    c:=x->-(evalf(Pi))*((lambda1*(1-x^2))+4*eta*lambda2*(1-x^3)+(lambda1*eta/5)*(1-x^5));
    rho:=x->rho0*int(c(x-x1)*rho(x1),x1=0..20)+rho0+((rho0)^2)*c(0);
    rho(x1):=rho0;
    for x from 0 by 1 to 20 do
      for i from 1 by 1 while shart=1 do
        t[i]:=rho(x)-rho(x1);
        if ((abs(sum('t[i]','i'=1..i)))^2<10^(-4)) then shart=1 else rho(x1):=rho(x) end if;
      end do;
      result(x):=rho(x);
    end do;
    with(plots):
    points:= { seq([x,result(x)],x=0..20) }:
    plot(points);
    and the mtlab code is:

    Code (Text):
    clear, clc, format long
    %
    n = 1000;
    a = 0;
    b = 20;
    h = (b - a)/(2*n);
    x = a:h:b;  
    %
    p0 = 0.57;
    E = (pi*p0)/6;
    L1 = (1 + 2*E)^2/(1 - E)^4;
    L2 = -((1 +.5* E)^2/(1 - E)^4);
    c0 = -1*pi*(L1+ 4*E*L2 +(E*L1)/5);
    %
    f = zeros(1,2*n + 1);
    p = zeros(1,2*n + 1);
    pp = zeros(1,2*n + 1);
    c = zeros(1,2*n + 1);
    p(:) = p0;
    e = 0.01;
    for i=1:100            % counts iteration steps
       for l = 1:2*n+1      % counts x
            for j=1:2*n+1    %counts x'
                c(j) = -1*pi*(L1*(1 - (x(l) - x(j))^2) ...
                    + 4*E*L2*(1 - (x(l) - x(j))^3) + E*(L1/5)*(1 - (x(l) - x(j))^5));
                f(j) = c(j)*p(j);
            end
            for k=3:2:2*n-1
                T1 =  sum(f(k));
            end
            for k=2:2:2*n
                T2 = sum(f(k));
            end
            m = h*(4*T1 + 2*T2 + f(1) + f(2*n+1))/3;
           pp(l) = p0*m + p0 - ((p0)^2)*c0;
       end
       for l = 1:2*n+1
           if abs(p(l)-pp(l)) > e
               break;
           end
       end
       if l == 2*n+1
           break;
       end
       p = pp;
    end
       for l = 1:2*n+1
           if abs(p(l)-pp(l)) > e
               break;
           end
       end
       if l == 2*n+1
           break;
       end
       p = pp;
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you help with the solution or looking for help too?
Draft saved Draft deleted



Similar Discussions: Problem with matlab or maple
Loading...