Problem with matlab or maple

  • Maple
  • Thread starter Majid
  • Start date
  • #1
23
0

Main Question or Discussion Point

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:
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:
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;
 

Answers and Replies

Related Threads on Problem with matlab or maple

Replies
4
Views
3K
Replies
2
Views
4K
Replies
2
Views
41K
Replies
9
Views
47K
Replies
3
Views
3K
Replies
2
Views
6K
Replies
17
Views
45K
Replies
18
Views
138K
Replies
10
Views
20K
Replies
2
Views
3K
Top